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ABSTRACT 

MINIMAL NORM CONSTRAINED INTERPOLATION 

Larry Dean Irvine 
Old Dominion University, 1985 
Director: Dr. Philip W. Smith 

In computational fluid dynamics and in CAD/CAM a physical boundary, 
usually known only discreetly (say, from a set of measurements), must 
often be approximated. An acceptable approximation must, of course, 
preserve the salient features of the data (convexity, concavity, etc.) 

In this dissertation we compute a smooth interpolant which is locally 
convex where the data are locally convex and is locally concave where 
the data are locally concave. 

Such an interpolant is found by posing and solving a minimization 
problem. The solution is a piecewise cubic polynomial. We actually 
solve this problem indirectly by using the Peano kernel theorem to 
recast this problem into an equivalent minimization problem having the 
second derivative of the interpolant as the solution. 

We are then led to solve a nonlinear system of equations. We 
show that with Newton's method we have an exceptionally attractive and 
efficient method for solving this nonlinear system of equations. 

We display examples of such interpol ants as well as convergence 
results obtained by using Newton's method. We list a FORTRAN program 
to compute these shape-preserving interpol ants. 

Next we consider the problem of computing the interpolant of 
minimal norm from a convex cone in a normed dual space. This is an 
extension of de Boor's work on minimal norm unconstrained interpolation. 



1 . The Natural Spline Interpolant 


We consider the problem of computing an interpolant to given data. 
Throughout our discussion we shall denote the data 

(t ,y ) 1 = 1,2, . . . ,n 

where a = t^ < < . . . < t = b and in this chapter we place no 

restrictions on the numbers y . There are, of course, many such 
mterpolants which we can form. For example, we can calculate the 
unique polynomial p of order n (degree n-1 or less) which interpolates 
the data. However, as pointed out in [deB(l), chapter 2], for large 
n (and especially for equally spaced points t ) the polynomial inter- 
polant is notorious for large changes m its first derivative near the 
endpoints. Figure (1.1) displays the polynomial interpolant to the 
function 

f(t) = 1 - -- s - ] 2 n - (7 77 t} 

at the points t^ = (i-l)/10 for i = 1,2,..., 11. Since 0 < y s 1 for 

each i, we expect a good interpolant to remain reasonably close to 

these bounds. However, because of its behavior near the endpoints, 

the polynomial interpolant fails to model the data well. This behavior 

is typical of high-order polynomial mterpolants. 

In order to decrease the unnaturally large changes in the first 

derivative characteristic of the polynomial interpolant, we wish to 

calculate the interpolant which "bends" the least over all suitable 

mterpolants. The norm of the second derivative of an interpolant 

will furnish a measure of the bending of the interpolant so we pose a 

( 2 ) 

minimization problem on [a,b], the Sobolev space of functions with 
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Figure (1.1): The Polynomial Interpolant. 
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second derivatives in the normed linear space I-^la.b]. Let A denote 
the set of all mterpoiants m the boboiev space. We consider the 
minimization problem 

Find f .,.0 A such that ||f # ^) || ^ < |) f ^ ^ ^ ^ for all ft A. (1.1) 

We shall see that the solution to (1.1) is piecewise cubic with two 
continuous derivatives; that is 


f*(t) = P i (t) if t i £ t < t i+1 

for l = l,2,...,n-l where p is a cubic polynomial and f # is m 

2 

C [a,b]. We follow the pattern m [deB(l), chapter 5], taking advan- 
tage of the fact that L 2 [a,b] is not only a normed linear space, but 
also a Hilbert space with an inner product defined by 


(f.g) = 


f(t)g(t)dt 


for any two elements f and g in L 2 [a,b]. 

Assume f is an element of A. (The set A is nonempty since it 

contains the polynomial interpolant . ) We shall use the Peano kernel 

( 2 ) 


theorem to obtain a set of equations for f 
Theorem of Calculus we have 


By the Fundamental 


tt 


f(t) = f (a) + 


f (1) (s)ds 


( 1 . 2 ) 


We 


ft (i) f 

integrate j f (s)ds by parts noting that udv = uv - Jvdu. 


Let 


u(s) = f^^(s) and dv(s) = di 


du(s) = f^^(s)ds and v(s) 


-(t-s) 


so that 
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where t is a constant. Hence 


r*u), 


:(D, 


'(s)ds = (t-a)f v '(a) + 
J a 

and so (1.2) becomes 


(t-s)f^ 2 \s)ds 


f(t) = q x (t) + 


- ( 1 ) j 


(t-s)f ^ 2 \s)ds 


(1.3) 


where q^(t) = f(a) + f' 1 '(a)(t-a). (This is actually a Taylor's series 
with integral remainder.) 

To acquire constant limit of integration we can write (1.3) as 

f(t) - q x (t) + f(t-s) + f (2)(s)ds (1.4) 

J a 

where (h) + , the positive part of the function h, is defined by 


(h) + (t) = 


h(t) if h(t) 6 0 

0 if h(t) S 0. 

Now we consider the divided difference operator. Given a function 


g and a set of points {t ,t t }, the m-th divided difference 

of g - denoted by [x^ . ,T i+m ]g( • ) ~ 1S the leading coefficient 

of the polynomial of order ra+1 which interpolates g at T i ,t ^ > • • * T 1+m 

(and hence is a function of T ,T ). The recursive relations 

i i+l i+m 


[ T p ] g( * ) = g(T p ) 


f T l ,T l+l’ 


»T 


i+m]g( * ) = 


[t 


1 + 1 


JLL 


T n-ml8 ( ' ) -i T 1 ..... T 1+ „-l)g ( - ) (1.5) 


T , - T 

i+m l 


hold if T ^ x (which we assume for our data). Presently we are 
interested in the case m=2. Equation (1.5) becomes (with T 1 = t i ) 
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(i 


i+2 




8(t i+2 ) ~ g(t i+l ) 
t I+2 ~ t i+l 


8(t i+l ) ~g (t i ) 
t x + l ~ C i 


( 1 . 6 ) 


which is computable for i=l,2, . . . ,n-2. 

Notice that [x^ ,T x+ ^, . . . , T ]p(*) = 0 if p is a polynomial of 
order m or less (degree m-1 or less). (From equation (1.6) we see 
that (t i+ 2 ~ t x )[ t^ , t i+ i , t i+ 2 ]g(* ) measures a difference m slopes; the 
difference in slopes being zero if g is linear.) 

Now we apply the (scaled) second-divided difference operator 
(t ^ + 2 ~ t 1 » t 1 +i » t 1 +2^ t0 (1*^) and interchange the order of the 
integral and divided difference operators to obtain 


where 



•b 

g(s)N (s)ds 
a 


1 = 1 , 2 , . . . ,n-2 


d i,2 (t i+2 t i^ t i ,t i+l ,t i+2^ f( ’ ) 


y i+2 - y i+l y i+l ~ y i 

t i+2 " t i+l t i+l ~ 


(1.7) 


( 1 . 8 ) 


N i, 2 ( ° = (t i,2 ~ t i )[t i* t i+r t i+2 ]( *- s) + 


(t „ - s) - (t . - s) (t - s) - (t - s) , 

v i +2 ' + v 1+1 ' + v 1+1 + 1 + 


fc i+2 t i+l 


t . - t 
l+l i 


(1.9) 

( 2 ) 

and g = f v . We call N „ the (normalized) linear B-spline (or 

l » ^ 

B-splme of order 2) with knots t , t . and t The graph of N „ 

1 1 + 1 1 + 2 or 1,2 


is displayed m figure (1.2). 
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Figure (1.2): The Normalized Linear B-splme. 
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We have shown that if f is an interpolant in the Sobolev space 

( 2 ) 

(f e A), then g = f v ' satisfies (1.7). Let the set B consist of all 
functions which are m L^Lajb] and which satisfy (1.7). 

Now consider the problem 

Find g., ; . e B such that || g. ;s .|| < jj g|| for all g Z B (1.10) 

A unique solution exists since (1.10) is a minimal norm problem over 

a nonempty closed convex set m a Hilbert space. Furthermore, the 

( 2 ) 

solutions of problems (1.1) and (1.10) are related via g^ = f* . 
Hence, to compute f.,,_ we can first calculate g^ and then integrate g^ 
twice. Since much of our emphasis will be on g^, rather than f^, we 
shall call g^. the interpolant of minimal norm. 

For brevity we denote the index m = n-2, the B-spline N = N „, 

X X 9 z. 

and the divided difference d = d „. We also define the vector-valued 

i i,2 

function TrL^Ca.b] -* R m by 

( b 

(Tx)^ = x(t)N i (t)dt i=l,2, ...,m. 

■'a 

To solve problem (1.10) we shall show that g^., the interpolant of 
minimal norm, is the intersection of two specific sets — one an ortho- 
gonal complement of a subspace and the other a translate of a subspace 
— in L 2 [a,b] via a variation of the Projection Theorem. If W is a 
closed subspace of a Hilbert space H and if x is an arbitrary element 
of H, then the Projection Theorem states that there exists a unique 
element w q in W satisfying 

|| x - w || S || x — w |j for all w e W (1.11) 


and characterized by 
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(x - w q ,w) = 0 for all w e W. 

Hence x - is in W "*■, the orthogonal complement of W. The proof of 
the Projection Theorem can be found m any book dealing with Hilbert 
spaces (for example, [L, page 517]). The next proposition will serve 
as the actual form of the Projection theorem which we shall use. 

Proposition ([L, page 64]): Let W be a closed subspace in a Hilbert 

space H. For a fixed element x _in_ H define V : = x + W. Then there 
exists a unique element x q ajl ^ of minimal norm. Furthermore , x q i_s 
m W 

(The translate V is called an affine set or linear variety.) Notice 

that x q is the intersection of the orthogonal complement of W and the 

translate V of W. In fact, (1.11) reveals that x = x - w . 

o o 

Define 


W: = {zcL^a.b] : Tz = 0} 

which is a closed subspace m L^tajb]. Let geL 2 [a,b] be any element 
such that Tg = 4* (Equivalently, let g be any element of B.) Then 
B = g + V/ and B corresponds to the linear variety m the proposition. 
Hence is the unique element m W 1 satisfying Tg^. = jd. 

We consider the contents of W 1 . Any element which is orthogonal 
to each is also orthogonal to any linear combination of the B- 
splmes. Hence S: = spa^N^,^, . . . ,N m ) is a subset of W 1 . We now 
show that W 1 is a subset of S (and hence S = W ’*’) by contradiction. 
Assume that there exists an element y which is in W but not in S. 
Since S is a closed subspace there exists (by the Projection Theorem) 
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an element s in S such that 
o 

j | y — s || S II y _ s || for all s e S 

with y - in the orthogonal complement of S. This implies 

T(y -s) = 0or(y-s)eW. However y - s is also m W 1 since 
o w o 3 o 

both y and s are in W 1 . Therefore (y - s ) =0 and S = W 1 . 

1 o J o 

In summary, g„. is characterized by 


m 

8* = Z « N 

i=l 

(since g # is m the span of the B-splmes) where the coefficients 

a. , a„, . . . ,a are chosen to satisfy 
l z m 

m 

( Z a N ,N ) = d i=l,2, ...,m (1.12) 

J=1 J J i i 


(since Tg^ = d^). Equation (1.13), a system of m linear equations m m 
unknowns, can be written m matrix notation as 


A a = d 


(1.13) 


where the symmetric matrix A has entries A = (N^.N^). 

Because the B-splines are linearly independent, the matrix A, a 
Grahm matrix, is nonsingular and hence a unique solution exists for 
any given _d. Furthermore, since N has support [t ,t _], the matrix 

1 X x+z 

A is tridiagonal. For any _xeR m we have 
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x Ax = Z x (Ax) 
“ “ i=l 1 “ 1 


m m 

Z x (N , Z x N ) 
i=l 1 1 j«l 2 2 


m m 

( Z x N , Z x N ) 
i=l 1 1 j=l 2 2 

m 


Z x N _ 

l l " 2 

1=1 


S 0 


with equality holding if and only if jc = 0. The matrix A is hence 
positive definite and (1.13) can be solved by Gauss elimination with- 
out pivoting, or, better still, by Cholesky decomposition. 

We note also that 

II g* II = a T Aa = a T d. 

The entry A , the integral of the product of two piecewise linear 
polynomials, can be computed exactly by Simpson's rule applied on each 
subinterval [t^,t^ + ^]. Denoting At^: = t^ + ^ - t^ and the midpoint 
of the interval [t^,t^ + ^] we have for i=l,2,...,m 


li 


'l+l 


N (t) dt + 
i 


'i+2 


N (t) dt 

i 


'i+l 


= (At i+1 /6)[N i (t i ) 2 + 4N i (z i ) 2 + N i (t i+1 ) 2 ] 


+ (At /6)[N (t ) 2 + 4N (z ) 2 + N (t ) 2 ] 

i+2 i i+l i i+l i i+2 


■ (t i+2 ~ t i )/3 ' 
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We also compute for i=l ,2 , . . . ,m-l 


A . = A . 
i,i+l i+l,i 


N iNi+i(t)dt 


(t i+2 t i+l^ /6. 


The solution g.„., being a linear combination of linear B-splmes, 
is piecewise linear (and continuous) with knots t . After integrating 
g.^ twice and applying the interpolation conditions, we obtain f # which 
is piecewise cubic (with knots t ) with two continuous derivatives. 

Define 3 £ R n via 


3 1 = < Vi 


i-2 , 3, . . . ,n 1 


and A3 = 3 1+ ^ - 3 . On [t ,t f* is defined by a unique cubic 
polynomial p.„ and hence f.„. can be determined by specifying the numbers 
p #i ^(t ) for i=l , 2 , . . . , n-1 and j=0,l,2,3. Then 


f*(t) = 


P*1 p*^) 


p., { t )(t-t ) + p,. ; (t )(t-t )' 

+ 1 1 1 1 1 1 


for t£ [t.^ t i+ ^]. Of course, (1.14) can be more efficiently evaluated 


by using nested multiplication. 
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The polynomial p # solves the differential equation 


P #1 (2) (t) = + (Ae i /At i )(t-t i ) 


(1.15) 


on the interval [t^.t^^] with boundary conditions p # (t^) = y^ and 
P *i (t i+1 ) = y i+l‘ Therefore 






1 


+ c (t-t ) + e 
i l 


i 


(1.16) 


where the constants c and e are evaluated as e 

ii l 


y and 
J i 


Ay. 


c = 


i At 




2 6 


At 


with Ay = y i+1 - y . From (1.17) we obtain 

p *i ( 0 )(t i ) ■ y i 


(1.17) 


P* ^(t ) = c 

t -^1 1 


(1.18) 


■s ( 2 V- s x 




where c is given by (1.17). A complete FORTRAN program for computing 

the natural cubic spline interpolant is given in Appendix A. 

Figure (1.3) displays the natural cubic spline interpolant that 

is m contrast to the polynomial interpolant of figure (1.1). 

We complete this chapter by posing (and solving) a generalization 

of problem (1.1). For k fixed satisfying 2 < k < n, let A(k) be the 

(k) 

be the set of interpolants (to the data) which are in L), [a,b]. We 
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Figure (1.3): The Natural Cubic Spline Intemolant. 
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consider the problem 


Find £ A(k) such that || 2 = Hf^^ i| 2 f° r all 


Let f be an element of A(k). Since (1.3) is valid for 
integrate by parts again (assuming k >2) to obtain 

• t 

I 

f(t) = q 2 (t) + 


(t ~ S> f (3) (s)ds 


2' 


J a 


where 


q 2 (t) = f (a) + f^\a)(t-a) + f (t-a) 2 . 


In general, after integrating by parts k-1 times we obtain 


«t) - 


Vi (t > + 


k-1 

(t-s) k %(k) 
(k-1)' 1 


(s)ds 


or 


f(t) = 


Vi (t) + 


V ^~1 

(t ~ s) + f (k), . , 

Wiyr f (s)ds - 


Now we take the (scaled) k-th divided difference of (1. 
obtain 

fb 


i,k 


g(s)Ni k ( s )ds i=l,2, ... ,n-l 


f e A(k) 

(1.19) 

f, we can 

( 1 . 20 ) 


( 1 . 21 ) 


( 1 . 22 ) 

22) to 

(1.23) 


a 
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where 


d i,k ■ (k - 1),(t i+k- t i )[t x, 


t ]f ( 
, i+k J v 


(1.24) 


N i,k< 8 > " 


(t . -t )[ t 
l+k i /l i, 


,w<- s > + 


k-1 


(1.25) 


( 2 ) 

(the normalized B-splxne of order k), and g = f 

Let B(k) denote the set of elements (m L 2 [a,b]) wich satisfy 
(1.23). Then the solution f^. to (1.20) is related to the solution 
to the problem 


Find g_ ;5 . e B(k) such that !|g. ;: / k ^|| 9 ^ ||g^ II 2 for § eB ( k ) (1*26) 


f ^ 

via = f.,,. ' . Furthermore, for some a e R we have 


n-k 

g„ = Z a N , . 

' j-i j J' k 


The coefficients a 1t a o ,...,0 , are chosen to solve the linear system 

1 2 n-k 

of n-k equations in n-k unknowns represented by the matrix equation 
A a = d_ where A is symmetric and positive definite with entries 


A 

ij 




Since g. ;s . is a linear combination of piecewise polynomials of order 
k, f^._ will be a piecewise polynomial of order 2k. We call f* the 
natural spline interpolant of order 2k. 



2 . A Minimal Norm Interpolation Problem 


in the Lp [a,b] Spaces 
For p such that 1 < p 5 °° we define the set 


: *[ 8 


G(p) : = <g £ L [a,b]: 


g(t)<t i (t)dt = 


'a 

for i=l 


» 2 , . . . , n^j 


8 o(t) < J,i(t) dt 


( 2 . 1 ) 


where { is a set of elements m L [a,b], q is conjugate to p 
(p+q = pq if p*°° and q=l if p= °°), and g Q is a fixed element of 
L [a,b]. Consider the problem 


Find g* £ G(p) such that ||g*|| p S || g || for all g e G(p). (2.2) 


In chapter 1 we solved (2.2) for the special case p=2; finding 
from a linear variety in a Hilbert space the element of minimal norm. 
The Projection Theorem came m handy to characterize g* as well as to 
guarantee uniqueness. However, for p *2 L p [a,b] does not have the 
orthogonality properties of a Hilbert space and hence, we cannot use 
the Projection Theorem to solve (2.2). Instead we solve (2.2) in this 
chapter by utilizing the Hahn-Banach theorem to characterize g. ;s . 

Uniqueness follows in the case 1 < p < °° by the strict convexity of the 
norm. This chapter, modeled after [deB(2)], motivates the use of the 
Hahn-Banach theorem in chapter 5. 

Let A be the linear functional defined on the subspace 


S: = span((t>j , . . . , <)) ) 


16 
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via 


A( Z , a <J> ) 
i=l 11 


(i=i 01 A)( t ) 8 o( t ) dt ' 


(2.3) 


Any element of G(p) (including g Q ) will serve as an extension of A to 
a bounded linear functional defined on all of L^[a,b]. Hence, 


A !l s s II 8 || p for all geG(p). 


(2.4) 


Conversely, any extension of A to a bounded linear functional defined 
on all of L [a,b], being identical to A on S, is represented by an 
element of G(p). 

The Hahn-Banach theorem guarantees the existence of an element 
g £ G(p) such that 


f(t)g(t)dt S || A || * |j f I! for all f eL [a,b] 
s q q 


This implies that jjg|| S ||A|| s which, taken along with (2.4), gives us 
|| g 'll = || A || and, therefore, a solution to (2.2). Now we characterize 


Let E a (j) be an element such that 
i=l 1 1 


II £ II = 1 and A ( z = ll A IL‘ 

1 =] 4 1=1 


(This element is unique if 1 < p < w since the norm is strictly 
convex.) Then 
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s 


n 

= X( £ a <t> ) 

i=l 1 1 


b 


n 


J a 


( 2 a <J> )(t)g(t)dt 
i=l 


S II l a * INI g || 

n=i y 


= llgl 


Therefore, equality holds throughout and we have 


n 

£ 

i=l 


( £ a i (D i )(t)g(t)dt = (I Z cmDJI q *||g ll p - 


A 

Since g and £ a cp are aligned, we must have 
i=l 


g(t) = || A|| • || £ a d> || q 1 signum ( £ a 4 )(t) 


i=l 


i i 


i=l 


i i 


We close this chapter by stating the interpolation problem that 

goes along with solving (2.2). Let p be a number such that 1 < p < 00 , 

(k) 

let k be an integer such that k 2 2, and let f e L v / [a,bl. Define 

op 

the sets 


F: = {ft Lp^k)[ a ,b]: f(t^) = f^Ct^) i=l,2,...,n) 
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and 


G: = {g 


e L [a, b] : 


g(t) N i k (t)dt = d 


x,k 


i=l,2, 


Then the problems 


Find f. ;; . £ F such that 



5 ||f (k) || p for all 


and 


Find f # e G such that || g^, II p = || g || for all g £ 


are equivalent and 


g*(t) = f* (k) (t) = 


n-k 

1 3 i N i k 
i=l 1 1,k 


q-1 n-k 

signum ( E BN 

i=l 1 1,R 


• , n-k} 


f e F 


G 


)(t). 



3 . The Convex Spline Interpolant 


The data { ( t ,y )} n _, are called convex if the point (t ,y ) 

1 1 1-1 1 rt 1 Cl 


lies on or beneath the line joining the points (t ,y ) and 

X 1 X 1 

(t ,y ) whenever 1 S l.. < i„ < i„ < n. Equivalently, 

1 3 1 3 ± Z J 


[t ,t ,t )f(*) = 0 
X 1 2 3 


(where f is any interpolant to the data) or 


d 

l 


y i+2 - y i + l y i+l ~ y i : 0 

fc i+2 ^l+l t i+l t i 


for i = l,2,...,m(= n-2). 

In this chapter we address the problem of finding, for convex 
data, the smoothest convex interpolant; that is, the convex interpolant 
having second derivative of minimal norm over all smooth convex mter- 
polants. The natural cubic spline interpolant, the smoothest of all 
interpolants, regrettably does not always preserve the convexity of 
the data. In chapter 1 we showed that f #t the natural cubic spline 
interpolant, has second derivative 


f * 


( 2 ) 


m 

E a N 


1=1 


J J 


where the coefficients • • • i 0 ^ satisfy (1.13). If any a x is 

negative, then f. ;s . is actually concave on a subset of [a,b]. 

Let {(t jy^)}^^ denote convex data and let A denote the set 

( 2 ) 

of convex interpolants m [a,b]. We assume that A is nonempty. 


20 
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(There are convex data for which A is empty. For example, let 

y = | t | and t. = -2, t„ = -1, t„ = 0, t, = 1, and t c = 2. The only 
i i 1 2 3 4 5 

/ o \ 

convex mterpolant is f(x) = |x|, which is not in [-2,2].) 

Using the Peano kernel theorem as we did in chapter 1 we can show 
that if f e A then T(f^)) = jf where TrL^Ca.b] ->R m is given by 
( T g 1 ) : = (g.N^) . Hence if 

B = {geL 2 [a,b]: g S 0 and T g = dj, 

then problems 

Find f # e A such that ||f w ^ || 2 £ II II 2 f° r f e ^ (3.1) 

and 

Find g.». E B such that ||g.^|| 9 £ ||g || 2 for all g £ B (3.2) 

(2) 

are equivalent and the solutions are related via = f# . Since B 
is a nonempty closed convex set, we consider (3.2) as finding the 
distance from a point to a closed convex set m a Hilbert space. 

Proposition ([L, page 69]): Let x be an element of a Hilbert space H 

and let K be a nonempty closed convex subset of H. Then there exists 
a unique element k q £ K such that 

|| x - k || S || x - k 1 1 for all keK 

Furthermore , k^ is characterized by 

(x-k.k-k)SO for all keK. 
o o 

Since we wish to find the element of minimal norm in B, X corre- 
sponds to the zero function and hence g^ is characterized by 


(g*>g-g-;.-) = 0 for a11 g e B. 


(3.3) 
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Propostion ([MSSW, proposition 2.1]): If there exist coefficients 

a^.o^, • • • ,a m satisfying 


b 
• a 


m 


( Z a N ) N (t)dt = d 
J=1 J J + i i 


i=l,2. 


,m. 


(3.4) 


then 


g fl = ( I a N ) . Furthermore, such coefficients exist i 
_ J=1 J J m 


f there 


exists g e B such that {I'M ^ are linearly independent over the support 
of g. 


Proof : Assume . . . ,ot m satisfy (3.4). 

assume g£ B, Define (h)_ = (-h) + so that 

h= (h) + - (h)_. 


Denote s = 


m 

Z a N 


J=1 


J J 


and 


Then 

((s) + , g-(s) + ) 

= (s + (s)_, g - (s) + ) 


= (s,g - (s) + ) + ((s)_,g - (s) + ) 
= ((s)_,g) - ((s)_, (s) + ) 

= ((s)_,g) 


a o. 

The last inequality is valid since both (s)_ and g are nonnegative 

functions. Hence (s) + satisfies (3.3). 

We now show that we can find coefficients a, ,a nf . . . .a so that 

12 m 

(3.4) holds by following the procedure employed in [MSSW]. 
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We begin by considering the problem 


inf 


m 


J a 


(Za N ) (t)dt : Zad =1 
j,U J + 


(3.5) 


and showing that if the infirnum is attained at some a, then for some 
positive constant C the coefficients Ca^, Ca m satisfy (3.4). 

If the mfimum of (3.5) is attained at a, then oi is a critical 
point of the Largrangian 


rb 


L(a , A) = 


m m 

( Z a N ) (t)dt + A(l- Zad) 
j=l J J + j-l 3 3 


(3.6) 


At a minimum of L we must have 


0 = 2 


J a 


( Z a N ) N (t)dt 
J=1 J J + i 


Ad i=l,2, 
i 


,m 


(3.7) 


and a • _d = 1 for some A. 

Now multiply (3.7) by a^ and sum over i=l,2,...,m to obtain 

rb 

mm m 

( Z aN ) ( Z aN )(t)dt - AZ ad =0 
1 J J + 1 i i ■, 1 1 

j=l J J i=l i=l 


or 


rb 


A = 2 


( Z a N X+ (t)dt 2 0. 
J=1 J J 


(3.8) 


If A>0, then (3.7) reveals that 
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m # 

( L a N ) N (t)dt = d 
! J J + i i 


J a 


1 = 1 , 2 , 


, m 


(3.9) 


where a = 2 a /A. 
J J 


If A = 0, then (3.8) reveals that 


fb 

m 

( Z a N ) , (t)dt = 0 
H JJ + 

a 

m 

where ot • ji = 1. This implies that ( Z a .N ) S 0. However, for any 

J=1 J J 

g e B we have 

^ m m 

( s a N )g(t)dt = l a (N ,g) 

j=l J J J= 1 J J 

a 


m 

= Z a d 


J=1 


J J 


= 1 

which is impossible because g is nonnegative on [a,b]. We conclude 
that A is strictly positive and, if the infimum m (3.5) is attained 
by some a, that (3.4) is solvable. We now show that the infimum is 
attained . 

Let {a^}^ be a minimizing sequence. If { ||a^^ || } ^_^ is 

2 

unbounded, then divide the objective function of (3.6) by |)a || and 
the constraint by ||a . There then exists a such that 
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| a || = 1 , 

1 Of) 


a • _d = 0, and 

■b m „ 

( X a N ) (t)dt = 0. 

j=l J J + 

J a J 

m 

We conclude that Z a N is nonpositive, but not identically 
J=1 J J 

/N 

zero. Since we have assumed there exists geB such that the 
B-splines are linearly independent on the support of g, 

m m 

0= lad = Z a (g,N ) 

J=1 J J j-1 J J 


m 

= (g, 2 a N ) 

J=1 3 3 

< 0 

which is a contradiction. Hence a minimizing sequence must be bounded 

and the infimum is attained via a convergent subsequence. This 

completes the proof of the proposition. 

We note that the existence of g e B, such that { N } m are linearly 

1 1=1 

independent over the support of g, in the previous proposition is 
guaranteed if d > 0 for each i. Then each g eB must be positive on 

some subinterval of t * ie su PP ort °f N i> for each l. 

Now we consider the implication of allowing d^ = 0 for some k. 

T 

As a specific example let t = (l-l) for i=l,2,3,4 and let = (1,0) . 

If g£. is the positive part of a linear combination of B-splines, then 
there must exist numbers and a ^ satisfying 
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3 

(a 1 N 1 + a 2 N 2 ) + N i (t)dt = 

•'0 


1 1 if i=l 
< 

( 0 if i=2 


(3.10) 


which implies that 0^ = - 00 . This is equivalent to the solution 
being identically zero on [1,3]. In fact, any g£ B must be of the 
form g = g X[q -jj- It is shown in [MSSW, theorem 3.1] that the 
solution to (3.2) is 

m 

g„ = ( £ a N ) x r 

J=1 

for appropriate coefficients a , . . . ,0^ where 


F: = [a,b]/{U (t ,t ) : d = 0} . 

J = 1 J J 

Hence the solution to (3.2) with t = (l-l) for i=l,2,3,4 and 
d = (1,0) T is 

= 3N 1 X [0,1]. 


Unless otherwise stated we assume d > 0 for each l for the remainder 

l 

of this chapter. 

Before we consider how Lo compute the coefficients 0t^, a . ..» a m 
which satisfy (3.4), we give a procedure for integrating g^. 

Define 3 X » A3 X , At^, an d Ay x as m chapter 1. We integrate g^«. on each 
subinterval separately, forming a piecewise polynomial, by 

solving the differential equation 

(2) ^i 

P*i (t > " +AF (t - t i» + 

1 

for t i 5 t » t ^ with boundary conditions p^Ct^) = and 


p * <t 1 +i ) * y i + i- 


(3.11) 
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Two integrations gives us 


n x At A3 0 

p*i (t) = id (6 i + at (t " t i )) + + c i 


(3.12) 


and 


P *i (t) 


At 

1 

6(A3 i ) 2 


(3 X + A3 i (t-t i )) 2 +c i (t-t i ) + e i 


(3.13) 


for constants c and e^. We proceed by cases. 

Case ] occurs when both 3 1 and 3 1+ ^ are nonnegative. The nonnega- 
tivity constraint is not active in this case and so (3.13) is equiva- 
lent to (1.16), although with modified constants c^ and e . The values 
p * 1 ^( t i) for J = 0,1, 2, 3, are given by (1.18). 

Case 2 occurs when 8 <0 and 3 . > 0. In this case p„ can be 

i l+l K *i 

defined by two polynomials: a linear polynomial q ^ defined on 

[t^jT ] - where the nonnegativity constraint is active and hence the 

second derivative is zero - and a cubic polynomial defined on 

[t , t . ] where 
L l l+l 


t = t - 3 At /A3 

i i iii 


Applying the boundary condition p. ;i . (t^) = y we obtain e 
Applying p # (t ^) = we get an equation for c : 


(At ) 


6(AB l ) 


r~T 2 ( 3 i+i } + c i At i + y x = y 1+r 


(3.14) 

= V 


Solving for c we have 
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^i (B 1+ 1 } At i 

Cl " At x' 2(A3 1 ) 2 


(3.15) 


From (3.11), (3.12), and (3.13) we obtain 

q n (t ) = y 

^il i •'i 


q^P(t ) = c 

n il l l 


q ii^ t i> ■ 0 


■ 0 


q „(t ) = c (t -t ) + y 
H i2 l i i i J i 


q ( ^(x ) = c 

n i2 x l 


(3.16) 


q i2 >( V ' 0 


q^)( Ti ) = A3 i /At i 

where T and are given by (3.14) and (3.15) respectively. 

Case 3 occurs when 3 >0 and 3 , < 0. In this case p„ is 

i l+l K *i 

defined by a cubic polynomial q on [t^.T^] and by a linear poly- 
nomial q ^2 on [ T 1> t i+ ^] with T defined by (3.14). These polynomials 
are determined by the values 
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= y i 


q^ )( ti } = c x 4- (6 i ) 2 At i /(2A3 i ) 




(3.17) 


q il )(t i } = A3 1 /At i 


q., 0 (x ) = c (t -t ) + e 

1 2 x ill i 


q ( i } (T ) = c 
^i2 i l 


= 0 


^2 >(t i) " 0 


where and e are given by 


Ay 1 (3 1 ) At x 

Cl At i 2(A3 i ) 2 


(3 l ) 3 (At i ) 2 
8l Yl 6(A3 x ) 2 


Case 4 occurs when 3 X and g ^ are both nonpositive. In this 
case we obtain a linear polynomial defined on ft^.t an ^ determined by 
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p.„ (t ) = y 
,r i i 1 




p *i >(t 1 ) * 0 


pif (t,) = 0. 


(3.18) 


Since g # is piecewise linear and continuous (with knots at the 
t 's and T^'s), f^ will be piecewise cubic with two continuous 
derivatives (if d > 0 for each l). We call f # the convex cubic 
spline mterpolant. 

Now we turn our attention to the task of numerically calculating 
the coefficients , 0 ( 2 > • • • »ct m which satisfy (3.4). We continue to 

rp 

assume that d >0 for each i. Define F:R ra -»- R m by F = (F , F„,...,F ) 
i ■'12m 

where 


F x (a) 


,b 

m 

( I a N ) N (t)dt 1=1,2, 
j=l ] 3 
■'a 


, m. 


(3.19) 


We wish to solve F(x) = d. 


One method is to use Jacobi iteration. An initial guess 


(o) 


X' ' = (x ' ' , X 


( 0 ) 


2 ^°),..., 1S chosen and a sequence 


is generated by calculating x^ k+ ^, once _x^ k ^ is known, by solving 


x (k) 
i vX l 


F ( x ; ' , . . . ,x; ;,x' ' ,x 


(k) (k+1) (k) 

i-l’ i ’ i+l’***’ 


x< k >) = d 

m l 
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f (k+1) r 

tor x for each 

involves calculating 


l. A modification 


(k+1) (k) 

x , once x 


the Gauss-Seidel iteration, 
is known, by solving 


F 

1 


(x< k+1> 


X 


(k+1) (k+1) (k) 

1-1 * 1 ’ 1 + 1 ’ 



= d 


l 


for x (k+1) for i=l, 2,..., m in succession. Both Jacobi and Gauss- 

l 

Siedel iterations converge globally as proved in [IMS] 

Now we consider Newton's method to solve G(x) = F(_x) - _d = 0 . 

We pick a suitable initial guess _x ' and form a sequence] x _ 
by solving 

(VG)(x (k) )(x (k+1) - x (k) ) = ~G( 2 i (k) ) (3.20) 

(k+1) (k) 

for x/ once x. is known. Since VG = VF, we can express (3.20) 
alternately as 

(VF)(x (k) )(x (k+1) ) - x (k) ) = d - F(x (k) ). (3.21) 


The entries of the Jabocian matrix VF are 

rb 


(VF) xj (a) = 


m 


( Z a,N.)°NN (t)dt 
k= i k k + J 1 


(3.22) 


m 


where ( I a N ) is the characteristic function for the support of 
k=i k k + 


m 

( Z a N ) + . We see that VF is symmetric and tndiagonal at each a. 
j=i 3 3 

We now characterize those a for which (VF)(a) is positive definite. 


Lemma (3.1) : The Jacobian (Vf)(a) is positive definite if and only if 
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( Z a, N, ) does not vanish identically on any of the subintervals 
k=l R R + 




Proof : 


For any x e K we have 


m 


_x (VF)(a)_x = I x 
i=l 1 


Z (VF) (a)x 
1=1 1J ~ J 


f b 

J a 


m mm 

( Z a N )°( Z x N )( Z x N )(t)dt 
k=l k k + j=l 1 3 i=l 1 1 


r b m m 

( Z ct k N k )°( Z x N r(t)dt 

k=l i=l 

a 


i 0 


,] for each l, 


If ( £ a N ) does not vanish identically on [t ,t „ 
j=l J J + 1 1+/ 

then equality holds if and only if x =0 for each l. If there exists 


some k such that ( Z a N ) is identically zero on [t, ,t, _], 

1 11“]" K Kt Z. 

j 


then 


equality does hold for the nonzero vector x. defined by x^ = 6^ k 

for each l. This completes the proof of the lemma. 

From (3.20) we see that 


F i _(o) 


m 

Z a 


1=1 


J 


■b m 

( k f x \Nk)>jNi(t) d t 

J a 


m 

= Z a (VF) (a) 
1=1 3 1J 
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so that F(u) = (VF)(a) a. Newton's method - equation (3.22) - takes 
the form 


(VF)(x (k) )x (k+1) = d. 


(3.23) 


Theorem (3.2) : If (vF)(x^^) is positive definite, then (VF)(x_^ <+ ^) 

is positive definite for each k and, hence, Newton's method - equation 


3.23) - is always well-defined. 


Proof : Having the known values x ' , we wish to determine the values 

(k+i) . . 

x satisfying 


S(k) J =1 J 


( ™x (k+1) N )N (t)dt = d i=l,2 m (3.24) 


where S(k) is the support of ( E x ) . Since (VF)(x5 k ^) is 

1=1 J J + 

positive definite, then S(k) U [ t^ , t^^ ] contains an interval for each i. 

Since d >0, then ( Z x ^ ) is positive on some subinterval of 

1 i J J + 

J=1 

[t^jt 2 !* Hence, (7 F) (x^^^ ) is positive definite. This completes 
the proof of the Theorem. 

Note that if has all positive components (for example, if 

(0) m (1) 
x = 1 for each 1 , then S(0) = [a,b] and Ex N is the second 
1 1 J J 

J=1 

derivative of the natural cubic spline mterpolant. 

Now we assume that d^ = 0 for some k. In this case special care 

must be exercised since {x _ n may diverge to -°°, preventing any 

k j - U 

numerical convergence. We already know that d^ = 0 implies that the 
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data points (t k >y k )» ^ t k+l ’ y k+l^ ’ and ^ L k+2 ,y k+2^ are colllnear and > 
hence, any convex interpolant must be linear on [ t ^ , t ^ 2 1 * ^H ulva_ 
lently, the second derivative of any convex interpolant must be zero 


OH U k .t k+ 2 J. 


Hence g.„ is of the form 


ViVA lx l"V 


' [t k+2' bl 


Since the value of x^ is immaterial and the k-th equation is 
automatically satisfied, the number of equations and unknowns each 
reduce by one. For computational convenience (3.23) can still be used 
with the following modifications: (VF)^ = 1, (VF) k = 0, and 

- »• 


If d^ = 0, then the solution is discontinuous at t^ if x^_-^ > 0 
and is discontinuous at t ^ + 2 1 f x k+i > 0. df the solution is discon- 
tinuous, then f.,,. will have only one continuous derivative. 

A further problem is encountered when d^_^ and are both 

zero, but d^ is nonzero for some k. Any nonnegative function g which 
satisfies the (k-l)-st and (k+l)-st equations can not satisfy the k-th 
equation since g is identically zero on [t^_^,t^ + ^] and on (t L , x1 »t^ x<J ] . 


•k+1* k+2- 


We conclude that there does not exist any convex interpolant in 
( 2 ) 

y [a,b] (and no solution to the problem as posed). However, we can 
find a convex interpolant whose second derivative is of the form 
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( Z X N ) { X r t 1 + X r , Ul) 
j=l J j 1 ,t k-l J lL k+3 ,bl 

satisfying all but the k-th equation. We already know that this convex 

interpolant must be linear on f t, ..t, ,] and on [t. , ,t, „] and, 

k-1 k+1 k+1 k+3 

hence, piecewise linear on [ t^_^ , t^ ^ ] . If d^ is nonzero, then there 
will be a discontinuity in slope at t^ For the convenience of 
utilizing (3.23) we can set d, to be zero to satisfy the k-th equation. 


The discontinuity in slope will show up after we integrate the solution 
to obtain the interpolant. 

Figure (3.1) displays the natural cubic spline interpolant to the 
function 


f(t) (0.05+t)(1.05-t) 

at the knots t-^ = 0, ^2 = 0.1. *-3 = 0.4, t^ = 0.7, tj. = 0.8, and 

tg = 1.0. Figure (3.2) displays the convex spline interpolant to this 

function. Table (3.1) shows the convergence results for Jacobi, Gauss- 
Seidel, and Newton's method iterations taken from [IMS]. Note the 
quadratic convergence characteristic of Newton's method. These conver- 
gence results are typical. 
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TABLE 3.1 
||F(x (n) ) - d || 2 


Iteration Gauss 


Number 

Jacobi 

Seidel 

Newton 

1 

.46 

X 

2 

10 

.27 

X 

io 2 

.19 

X 

io 2 

2 

.28 

X 

io 2 

.11 

X 

io 2 

.85 

X 

60 1 

3 

.75 

X 

10 1 

.42 

X 

io 1 

,29 

X 

io 1 

4 

.12 

X 

10 2 

.18 

X 

io 1 

.49 

X 

o 

o 

r— H 

5 

.26 

X 

io 1 

.75 

X 

10° 

.14 

X 

10 1 

6 

.49 

X 

io 1 

.31 

X 

10° 

.11 

X 

io ' 4 

7 

.10 

X 

io 1 

.13 

X 

10° 

.71 

X 

io ' 11 

8 

.21 

X 

io 1 

.55 

X 

io ' 1 

.49 

X 

10~ 12 

9 

.43 

X 

10° 

.23 

X 

1 — 1 
o 

1 

H-* 

- 

— 

— 

10 

.86 

X 

10° 

.96 

X 

io " 2 






20 

.11 

X 

10 _1 

.16 

X 

io -5 

- 

— 

— 

30 

.14 

X 

10 

.26 

X 

10~ 9 

- 

— 

— 




-5 



-13 




40 

.18 

X 

10 

.58 

X 

10 


— 

— 




-7 







50 

.24 

X 

10 

- 

— 

— 

- 

— 

— 

60 

.30 

X 

10“ 9 

- 


— 

- 

— 

— 




-11 







70 

.39 

X 

10 

- 


— 

- 

— 

— 



. The Shape-Preseryinc Spline Interpolant 


We addressed in chapter 3 the problem of finding, for convex data, 
the smoothest convex interpolant. We begin this chapter by considering 
the problem of finding, for concave data, the smoothest concave inter- 
polant. Then we continue the chapter by examining the problem of 
finding, for general data, the smoothest interpolant which is locally 
convex where the data are locally convex and is locally concave where 
the data are locally concave. 

Let { ( t ,y )} n , denote concave data and let A denote the set of 
1 i=l 

( 2 ) 

all concave interpolants in L ^ [a,b]. Assume A is nonempty. Using 

the Peano kernel theorem as we did m chapter 1, we see that, if f £ A, 


f (2) (t)N (t)dt = d 1=1,2, ... ,m(=n-2) 


Equivalently, we have T(f (2) )= d. 


Defining 


B: = {g e L„[a,b] : g i 0 and Tg = d_) , 


we conclude that the problems 


Find f # £ A such that ||f i ./ 2 ^|| ^ = ||f^|| 2 f° r all f £ A (4.1) 


(the problem of finding the smoothest concave interpolant) and 


Find g* £ B such that l[g. ; , || 9 < || g \\ ? for all g ' B 


38 
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( 2 ) 

are equivalent and the solutions are related via g^,. = f^. 

Of course, the smoothest concave interpolant to the concave data 
{(t^.y^)}^ is the negative of the smoothest convex interpolant to 
the convex data { ( t , -y )} . We highlight this with the following 

proposition . 


Proposition [MSSW]: If there exist coefficients . . .ct^ 

satisfying 

/■b m 

- ( £ a N )_N (t)dt = d i=l,2, ...,m (4.3) 

1 J J i i 

m 

then g { = - ( Z a N ) . Furthermore, such coefficients exist if there 
J= 1 J J - 

exists g € B such that ^ are linearly independent over the 

support of g . 


We note that the existence of g c B, such that { N^} are 

A 

linearly independent over the support of g, in the previous proposition 
is guaranteed if d < 0 for each 1 . Then each g e B is negative on 
some subinterval of [t^,t 2 ]* t ^ ie su PP ort for , for each l. 

Now we consider the problem of finding, for general data, a 
smooth shape-preserving interpolant - a smooth interpolant which is 
locally convex where the data are locally convex and is locally concave 
where the data are locally concave. Assuming for the moment that d^ 
is nonzero for each l, we define the sets 
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Tj : = {[t i ,t i+2 ] : d i > 0} , 


T 2 : = ; d i > 01 • 


°1 : ' T ]/ T 2 ’ 


s 2 : . T 2 /T, , 


and 


: = [a.bJ/C^U ^ 2 ) • 


Now we define the sets 


A: = {f e l 2 (2) [a,b] : f (2) X Q 2 0 , f (2) X fi £ 0, 


and f(t ) = y i=l,2,...,n} 


(which we assume is nonempty) and 


B: = {g e I- 2 [a,b] : gX^ 2 0, gx Q 2 0, and Tg = .d } . 

We conclude that the problems 

Find f.;,. £ A such that ||f^ " 2 ^ || 2 S || f^ 2 ^ || 2 for all f £ A (4.4) 


and 


Find g* £ B such that ||g # || 2 £ ||g || 2 for all g e B 


(4.5) 
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( 2 ) 

are equivalent, and g^ = l.,,. 

The following proposition gives the solution to (4.5). We see 
that has the character of the convex spline interpolant, 

has the character of the concave spline interpolant, and f*X^ has the 

character of the natural spline interpolant. 


Proposition [MSSW]: If there exists coefficents 


sati sf ying 


b 

m 


(< Vj N i )+X “ 

J = 1 J 


m 


( L a N ) 
j=i J J 



,a 


m 


m 

+ ( Z a N )x 0 1 N (t)dt = d i=l,2, . . . ,m 
j=l J J “3 1 1 

(4.6) 


then 



Q 


Furthermore, such coefficients exist if there exists g e B such that 
{ N } ™ , are linearly independent over the support of g. 


We note that the existence of g £ B, such that { N }™_j_ are 

/* 

linearly independent over the support of g, m the previous proposition 
is guaranteed if d is nonzero for each 1 . Then each g£ B is nonzero on 


r 
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some subinterval of [t x ,t 2 ]* the support of N , for each 1 . 

We now solve (4.6). Define F : R m -*■ R m where F = (F^.F^, . . . ,F m )^ 
is given 


F x (x) 


m 

( Z x 

O J =1 
Q 1 


N ) N (t)dt 


3 J'+ i 



m 

( E xN ) N (t)dt 
i J J ~ i 
3 = 1 


+ 



( Z x N )N (t)dt 
]=l jJ 1 


i=l , 2 , . . . ,m 


(4.7) 


We use Newton's method to solve F(a) = d_. 

, (0) , r (0) 

initial guess x. we produce a sequence l_x , 


Picking a suitable 
x^^ ,...,} by solving 


(VF)(x (k) )(x (k+1) - x (k ) = A ~ F(x (k) ) (4.8) 


for 


(k+1) 

X 


once x 


(k) 


is known. 


The Jacobian matrix has entries given 


by 


(?F) (a) 


f b 

P(«)N (t)N x (t)dt 
J a 


(4.9) 


where 
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m m 

P(a) - ( E.N )V + ( Ul )\ + X„ • (A.10) 

j=l J J 1 j=l J J 2 3 


From (4.9) we see thatVF is symmetric and tridiagonal at each . 
We also note that 


m m 

P(x)( Z x N ) = ( Z x N ) Xq 
j= 1 J J J= l J J + “i 



+ 


so that F(_x) = (VF)(x)2c and, hence, (4.8) reduces to 


(VF)(x (k) )x (k+1) = d. 


(4.11) 


The following lemma (with proof similar to its counterpart in 
chapter 3) characterizes those a for which (vF)(a) is positive definite. 


Lemma(4 , 1 ) : The Jacobian (VF)(a) is positive definite if and only if 


P(ot) does not vanish identically on any of the subintervals 


[ t 1 » t i+ 2 ] for i=l,2, . . ■ ,m . 
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The following theorem is modeled after theorem (3.2). 

Theorem (4.2) : If (VF)(_x^^) is, positive definite, then Newton's 

method - equation (4.10) - is always well-defined . 


Note that if x.^*"^ is given by x = signum (d x ) for each l, 
then P(_x^^) is the characteristic function for the interval [a,b] and 


m (1) 

Z x 'N is the second derivative of the natural cubic spline 

J-l J J 


interpolant. 

If d^ = 0 for some k, then we already know that any shape-preser- 
ving interpolant must be linear on [ t^ , t ^ + 2 ] • In fact any g £ B must 
satisfy 




The solution m this case is of the form 


8 * - h{X [a,t,J + X [t,_ 0 b] 1 


"k+2. 


m 

h = ( Z a N ) X„ 
j.i j j + n i 


m m 

- ( Z a N ) X 0 + ( Z a N )X n . 
].lO - a 2 j.l J J n 3 


where 



45 


Since the value of ot^ is immaterial - the k-Lh equation F^(ot) = of 
(4.11) being automatically satisfied - the number of equations and 
unknowns reduce by one each. For computational convenience we can 
still use (4.10) by setting (VF)^ = 1, (VF)^ = 0, and (VF)^ 

Once we solve F(a) = jd we proceed to integrate g.^ which is 
piecewise linear (but not necessarily continuous, even if d^ is non- 
zero for each k) to obtain f^ which is piecewise cubic. On the 
interval [t^t f* 1S given by the solution to the differential 
equation 

p (2) (t) = 6 + £8 /At )(t-t ) (4.12) 

*i l l l l 


for 


t s t 

l 


S fc i+l lf ^i^i+l^S* 


P 1 (2) (t) = ce x + (A6 i /At i )(t-t i )) + 


(4.13) 


for t < t < t . if [t ,t ,JC f!,, or 
l - - l+l l l+l 1’ 


P x (2) (t) = -($ x +(A3 i /At i )(t-t i )) 


(4.14) 


for t S t 5 t , if [t ,t . ] C with boundary conditions 
l l+l i l+l 2 


P (t ) = y and 


P i (t i+1 } = y 


l+l’ 
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The function is either a cubic polynomial or piecewise cubic 
given by two polynomials q^ and defined on separate subintervals 
of [t^.t ^]. The solution p^ to (4.11) is given by (1.18). The 

solution to (4.12) is, depending on signum (3^ and signum (8 i+1 ), 
given by (1.18), (3.16), (3.17), and (3.18). The solution to (4.13) 
is determined by (1.18) if B x 2 0 and 0 ^ 2 0, by (3.16) if 6 1 > 0 

and 3 1+1 < 0, by (3.17) if 0 < 0 and 0 x > 0, and by (3.18) if 

8 i = 0 and 0 i 0. 

Figures (4.1), (4.3), (4.5) and (4.7) display the natural cubic 
spline mterpolants to the given data. Figures (4.2), (4.4), (4.6), 
and (4.8) display the corresponding shape-preserving mterpolants. 
Tables (4.1), (4.2), (4.3), and (4.4) give convergence results for 
Newton's method. Note the quadratic convergence characteristic of 
Newton's method. 

Appendix B lists a FORTRAN program for computing the shape-preser- 
ving cubic spline mterpolant. 






Figure: (A. 3) : The Natural Cubic Spline Interpolant. 



Figure (4.4): The Shape-Preserving Cubic Spline Interpolant 




Figure (4.5): The Natural Cubic Spline Interpolant. 



Figure (4.6): The Shape-Preserving Cubic Spline Interpolant. 




'll 


Table 4.1 


Iteration Number 
1 
2 

3 

4 

5 

6 
7 


i|F(2c (n) ) -d|| 2 
0.13 x 10 1 
0.67 x 10° 
0.25 x 10° 
0.42 x 10 _1 
0.12 x 10 -2 
0.88 x 10~ 6 
0.58 x 10 -12 
0.64 x 10~ 13 


8 
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Tab] o 4.2 


Iteration Number 
1 
2 

3 

4 

5 

6 
7 


||F(x (n) ) - d || 2 

0.12 x 10 1 

0.56 x 10° 

0121 x 10° 

0.36 x 10 _1 

0.11 x 10“ 2 

0.85 x 10” 6 
-12 

0.54 x 10 
0.70 x 10 -13 


8 
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Table 4.3 


Iteration Number 
1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 


II F(_x (n) ) - d || 2 

0.24 x 10 1 
0.16 x 10 1 
0.12 x 10 1 
0.90 x 10° 
0.53 x 10° 
0.20 x 10° 
0.26 x lO -1 
0.42 x 10~ 3 
0.97 x Hf 7 
0.37 x 10 -12 
0.21 x 10 -12 
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Table 4.4 


Iteration Number 

1 

2 

3 

4 

5 

6 
7 


||F(x (n) - d || 2 

0.29 x 10 1 
0.13 x 10 1 
0.50 x 10° 
0.11 x 10° 
0.56 x 10~ 2 
0.16 x 10~ 4 
0.13 x 10~ 9 
0.26 x 10~ 12 


8 



Constrained M i n nn i /at i un in a Du.i I Space 


5. 

Let C bo a convex corio : ri a normed dual space X with predual Y. 
Assume are e l ements of Y and define T: X R n by 


Tx = (x(y 1 ),x(y 2 ) 


.,x(y n )) 


T 


Let B: = {x £ C : Tx = d} for a given vector d_. Consider the 
problem 


Find x. ;! , e B such that ||x. ;; .|| S ||x|| for all x e B (5.1) 

of which (1.10), (3.2), and (4.5) are special cases. In this chapter 
we study existence and characterization of solutions to (5.1). The 
follov'ing lemma gives sufficient conditions for existence of a solution. 

■ss* 

Lemma(5.1) : If B is nonempty, if C is weak closed, and if Y is 

separable, then there exists a solution to problem (5.1). 


Proof : Let y : = inf { 1 1 x | j : x e C and Tx = d } . Let { x^} be a 

sequence in C such that 


Tx = d 
n — 


(5.2) 


and 


II x n l! = Y + 1 / n 


(5.3) 
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for each n. Since Y is separable, by Alaoglu's theorem there exists 
a weak* convergent subsequence of { x^} with weak* limit x. Since C 
is weak* closed we have x £ C, from (5.2) we have Tx = ji, and from 
(5.3) we have ||x|| S y (and hence || x|| = y). This completes the 
proof of the lemma. 

Throughout this chapter we assume that B is nonempty, C is weak* 
closed, and Y is separable. Since x„. =9 if d = 9, we assume also 

that d_ * 9 . The following proposition gives us sufficient conditions 
for C being weak* closed. 

Proposition (5.2) : If C is normed closed and if Y is a reflexive 

space, then C is weak* closed. 

Proof : Assume { x n } is a sequence in C with weak* limit x. We want 
to show that x is m C. We do this by contradiction. If x is not an 
element of C, then there exists an element y (an element of both the 
dual and predual of X) which serves to separate x from C in the sense 
that 

x n (y)> K 

for each n and 

x(y) < K 

for some constant K. This implies that 

lim x n (y) * x(y) 

n oo 
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which is a contradiction. Therefore x e C and C is weak"* closed. 
This completes the proof of the proposition. 

For y > 0 we define the convex set G(y) C R n by 

G(y): = {Tx : x e C and || x || 2 y} . 

We now show that G(y) = YG(1) and G(y) is closed. 


Proposition (5.3) : For each y > 0 we have G(y) = YG(1). 


Proof : By definition 

G(y) = {Tx : x e C and ||x |j 2 y} 

= { Tx : ^ e C and || x/y || 2 1} 

= {T(x/y) : ~ e C and |jx/y|| 2 1} 

= y{Tw : w e C and ||w|| 2 1} 

= yG(i). 

Proposition (5.4) : The set G(l) is closed. 

Proof : Assume { z^} is a sequence m G(l) which converges to _z. We 

want to show that is an element of G(l). Equivalently, we want to 
show that x e C exists such that ||x || 2 1 and Tx = _z. 
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For each n there exists x £ C such that x s 1 and Tx = z . 

n n n — n 

By Alaoglu's theorem there exists a subsequence of || x^ || which 
converges weak* to some x £ C. Hence || x || 2 1 and Tx = This 
completes the proof of the proposition. 

We define 


Y : 


= inf( y : d^ £ G(y) } . 


(5.4) 


Equivalently , 


Y = inf{y : There exists x £ C such that 
Tx = d_ and || x |J 2 yl 

= inf { || x || : x £ C and Tx = d_} . (5.5) 

By lemma (5.1) we know that there exists x^. £ C such that 
|| x^ || = y and Tx_.,_ = jd. We call x # an interpolant of minima] norm. 

We now attempt to characterize x,_ via the Hahn-Banach theorem. 

We begin by defining a functional p : Y -*■ R by 

p(y) = sup { x(y) : x e C and || x || Si}. 

Notice that if C = X (the unconstrained problem), then p is the norm 
on Y. In general, since we are taking the supremum over a subset of 
the closed unit ball U m X, we have p(y) S ||y || for all y £ Y. 

Since 0 is an element of C, we have p2 0. In convex analysis p is 
called the support functional of the convex set (x £ C : ||x|| S 1}. 
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Since C is weak* closed, the supremum is attained at some element of 
{ x C C : |J x (| 2 1} ; that is, for any y e Y there exists an x (a 
function of y) such that x e C, |jx j| s 1, and p(y) = x(y). In fact 
we have ||x|j = 1 unless x = 0- The following two propositions reveal 
that p is continuous, subadditive, and positive homogeneous. 

Lemma (5. 5) : The functional p is continuous. 

Proof ; Assume y^ and y 2 are elements of Y and define y = y^ - y 2 - 
Let x be the element in {x e C : ||x || < 1} such that p(y 2 ) = x(y 2 ). 
Since |x(y)| s ||y||, we have 

x(y 2 ) - !|y|| S x(y 2 ) + x(y) 
or 

x(y 2 ) - !| y !! Sx( Yl ). 

Therefore , 

P (y 2 ) - || y IN P ( yi ). 

The elements y^ and y 2 can be interchanged to obtain 


p (y x ) - II y II = p (y 2 ) 
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and hence 


|p(y 1 ) - p(y 2 )l s II y x - y 2 II • 


Lemma (5.6) : The functional p is subadditive and positive homogeneous 


(hence convex) . 


Proof : Assume and y^ are in Y. To show that p is subadditive we 

must show that 


p (yjl + y 2 ) s p(y x ) + p(y 2 ) • 


By definition 


p (yf+y 2 ) = SU P fxCy-j+y-?) : x e C and ||x|| £ 1} 


£ sup (x(y^) : x E C and ||x|| S 1} 

+ sup (x(y 2 ) : x E c and ||x|| 5 l) 

= PCy-^) + p (y 2 ). 

Now assume a > 0 and y c Y. To show that p is positive homoge- 
neous we must show that 


p(ay) = ap(y). 
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By definition 

p(oty) = sup {x(ay) : x e C and ||x || S 1} 
= a • sup {x(y) : x e C and ||x|| S 1} 
= ap(y). 


This completes the proof of the lemma. 


As an example we compute p for the case C = {x e L [a,b]: x 2 0} 
where 1 < p < For an arbitrary element g in L^[a,b], the predual 
of Lp[a,b] where p + q = pq, we have for any f e C with ||f || S 1 by 
the Minkowski inequality 


,b 

f(t)g(t)dt 

.a 


fb 


< 


f (t)g + (t)dt 

.a 




q 


■+" q * 
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Assuming g + * 0, let 


f 


(g>r 1 / 



p' 


Then we have f £ C, || f || ^ = 1, and 


f b 

f(t)g(t)dt = ||g + || . 

a 


Hence 


p(g) = sup { f(t)g(t)dt : f £ C and ||f|| 5 1} 


3 + 11 q 


If g + = 0, then p(g) = 0. 


Lemma (5.7) : For all a £ R n we have 


n n 

E a d S y 'p( Lay). (5.6) 

-iii - i 

i=l i=l 


Proof : Since y = inf{ y : ji £ G(y)}, we have d_£G(y + e) for any 

£ > 0. Hence for every positive integer n there exists £ C 
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such that Tx = d and II x j| S y + 1/m. Therefore, for any a £ R n 

m — " m " 3 - 


n 


E a d = E a x (y ) 
,11 ,i m w i 

1=1 i=l 


n 

= x ( E a y ) 
m , li 
i=l 


s ||x m !| p( Jay} 

1 = 1 

n 

= ( Y + l/m)p( E ay), 

i=l 

Now let m-*- 00 to obtain (5.6). This completes the proof of the lemma. 

Since we know that G(Y*) is closed from proposition (5.4), we 
could have used x.,,. in place of m the proof of lemma (5.7). The 
next lemma states that there exists a nonzero vector S £ R n such that 
equality holds in (5.6) . 

Proposition (5.8) : There exists a vector J3 e R n such that ||3|| =1 

and 


}i 

3 • d = y" ( E By,). (5.7) 

~ 1=1 1 1 

Proof : The vector d_ is an element of G(y ) , but not an element of 

G(y-e) for any e > 0« Hence the closed convex set G(y - e) and the 
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vector can be strictly separated by a hyperplane. This implies 
the existence of a nonzero vector J3(e) such that 

6 (e) * X < 8(e) * _d 

for all _y_ e G(y - e) and without loss of generality we may assume 
that || 3(e) || = 1. Equivalently, we have 

6(e) • Tx < 6(e) • d 

and by the linearity of T 

n 

x( Z 3 (e)y ) < 8(e) • d 

i=l 1 1 

for all x e C such that ||x|| S y - e. Hence we obtain 

n 

(y" - e)p( £ 6 (e)y ) < 6(e) • d^ 

i=] 

We can take the limit as r + 0 to obtain a vector 6. such that ||J3 || = 1 

n 

Y 0( Z 3 y ) s 6 • d. 

1 = 1 11 “ ~ 


and 
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We have the reverse inequality from lemma (5.7) and therefore 


3 


A = Y p( £ 3 y ). 

1=1 


This completes the proof of the lemma. 


Let A be a linear functional defined on the subspace 


S : = span(y 1 ,y , . . . , y ) 


by 


n 

A( £ a 


n 


y ) = £ a d 

i 11 i 1 t 

i=l i=l 


so that (5.6) can now be written 


A(y) £ y p(y) for all y e S. 

The Hahn-Banach theorem scates that there exists an element 
m X such that 

w(y) = A(y) for all y e S 

and 

v. 

w(y) £ Y p (y) 


w 


(5.8) 


for all y e Y. 


(5.9) 



Theorem (5.9) : The Hahn-Banach extension w is an interpolant of 


minimal norm. 


Proof : From (5.8) we see that Tv/ = d_ so that w interpolates the data. 

To complete the proof we show that w e C and ||w|| = y . 

We show that w is in C by contradiction. Assume w is not an 

element of C. Since C is weak 41 ' closed, there exists an element y 

'o 

m Y which strictly separates w from C in the sense that 

w(y Q ) > x(y Q ) for all x £ C. (5.10) 

Since C is a cone we have Ax e C whenever A > 0 and x £ C. Hence 
(5.10) implies 


0 S x(y ) for all x £ C 
o 


(5.11) 


(or p(y ) = 0) and 


w(y Q ) > 0. (5.12) 

However, from (5.9) and (5.12) we have 

0 < w(y Q ) 2 y*p(y Q ) = 0 


which is a contradiction. Hence w must be an element of C. 
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Lastly, we show that |[w j| = y . We already know that 

y' C S || w || (5.13) 

since w £ B (w £ C and Tw = _d ) . Because p is bounded above by the 
norm on Y, (5.9) yields 


w(y) S y || y j| for all y £ Y 


and hence 


II w I| S y*. (5.14) 

-J4. 

Taken together, (5.13) and (5.14) imply that ||w|| = y . This completes 
the proof of the theorem. 


Recall that foi a given element y m Y there exists an element 
6 1 o 

x (a function of y ) in C such that p(y ) = x(y ). Furthermore, 
either ll x 0 ll = 1 or x q is the zero element. The following lemma 

will lead us to the conclusion that, if p is differentiable at y , 

then p'(y ) = x . 

3 o o 

Lemma (5.10) : Let f be a functional defined on a normed linear space 

Z. If f is differentiable at x e Z and if there exists a linear 
— 0 

functional X such that 
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f(z Q ) + A( z-z q ) = f(z) 


(5.15) 


for all z in some neighborhood of z q , then 


A = (Vf)(z o ). 


Proof : Lot z = z -t tu where t > 0 and u t Z. Inequality (5.15) 

yields 


f(z + tu) - f ( z ) 

A(u ) < 2 2 . . ( 5 . 3 . 6 ) 

Since (5.16) holds for all t > 0 (and sufficiently small) and for all 
u £ Z, we have AS(Vf)(z Q ). Substituting -u for u in (5.16) yields 


A (u) > 


f(z - tu) - f(z ) 
o o_ 

t 


(5.17) 


for all t > 0 (and sufficiently small) and for all u e Z. Taken 
together, (5.16) and (5.17) imply A = (Vf)(z Q ). 

Corollary (5.11) : If p is differentiable at y Q e Y, then p' (y ) = x q . 

Proof : Since p(Y q ) = x o (y Q ) an< ^ x Q (y) “ p(y) for all y £ Y, we have 

p + x o < - y ~ y 0 } s p( - y) 
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for all y £ Y. By the previous lemma we have p'(y ) = x q . This 
completes the proof of the corollary. 

Inequality (5.6) motivates the problem 


n 

inf{ p ( E a y ) : Of • d_ = 1 } . (5.18) 

a i=l 1 1 


Notice that if a is any vector satisfying a • d_ = 1 and if x is any 
element of B, then 


n n 

1 = £ a d = x( Z a y ) 

i=l 11 i=l 1 1 


5 If x || p( Z ay) 

i=l 


and hence 


p( Z a y ) 

i=l 11 



This implies that the mfimum is positive (and, in fact, is bounded 

'/<■ __ 2 _ 

below by (y ) . If the mfimum is attained at some a £ R and if p 

n 

is differential e at Z a y , then we are led to a solution to (5.1) as 

i=l 1 
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the next theorem reveals. 


Theorem (5.12) : If there exists a e R n such that a • d_ = 1 and 


n n 

p( Z a*y ) = inf { p( Z a y x ) : a • _d = 1} 
i=l 1 a i=l 1 


n 

■w 

and if p is differentiable at Z a y , then 

i-'i' 

i=l 


Y 


P' 


n 

( Z a 
i=l 




is an interpolant of minimal norm. 


Proof : Problem (5.18) has Lagrangian 


n n 

L(a,A) = p( Z^ot y ) - A( Z a d - 1). (5.19) 

!=1 


If there exists a solution a to (5.18), then there exists A so 
* * 

that (a ,A ) is a stationary point of (5.19). Hence 


x(y ) - A d = 0 

l l 


i=l,2, . . . ,n 


(5.20) 
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where x = p'( £ a y ), x £ C, ||x|| = 1, and ot * d_ = 1 ■ 


i=l 


i' i 


We first show that X > 0. Multiply (5.20) by and 


sum over i to obtain 


n „ .... n „ ... 

x( £ £y) = A £ = A • 

i=l l-l 


n ^ 

Since x = p'( £ ct^), we h a ve 
i=l 


n ... n .... 

x( l a y ) = p( £ a y ) 
, i J i . i i 

i=l i=l 


so that 


A = p(£oty)20 
i=l 1 1 
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Actually, we know that since the infimum is positive, we have A > 0. 
We can also show this by contradiction. If A =0, then 


n 

x( Z ot y ) 
, i-'i 
i=l 


S 0 for all x t C. 


(5.21) 


Let s be any interpolant in C. (We know that there exists an mter- 
polant m C since B is nonempty.) Then 


s( T a y ) = lad 
i 1 1 ,11 

i=l i=l 


1 


which contradicts 
Now we show that 
interpolant in C. 

Y S |j x || /A = 1/A 


(5.21). Therefore, A >0. 

■Jr Jr 

Ay =1. From (5.20) we see that x/A is an 
Hence 


or 


Y A S 1 (5.22) 


Let w be an interpolant of minimal norm satisfying (5.9). Then 
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w( E a y ) s Y p( £ a v )• 


i=l 


i i 


i = L 


i i 


Equivalently, we have 


n n 

w( L a" y ) 5 Y x( Z tt'y ) 
, 1 1 , 1-^3 

1=1 1=1 


which leads to 


y“* 


(5.23) 


Taken together, (5.22) and (5.23) imply 


1 



Tins concludes the proof of the theorem. 

We consider now the problem of determining when the mfimum is 
attained in (5.18). From proposition (5.8) we know that there exist 
a nonzero vector 8 such that 


-Yc n 

0s8*d=YP( Z B y ) 

i=l 11 
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If 3 *d_ > 0, then the infimum is attained in (5.18) at a = J3/(J3 * d.) . 


Proposition (5.13) : If _d is in the relative interior of 


S: = {_r : £ £ G(y) for some y } , 


then there exists a vector 3 such that 


y. ft 

1 = 3 • d = y' p( Z 3 y ). 

i=l 


Proof : We prove by contradiction. Assume that every vector 3. which 

satisfies 


■>£ ft 

3 • d = y p( I By) 

i=l 1 


also satisfies 3_ • d_ = 0. Without loss of generality it can be 
assumed that there exists a nonzero vector 3 such that 


% n 

0 = 8 • d = y p( Z 3 y ) 
i=l 1 1 


3 • £ 2 0 for all £ G G (y ). 


and 
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In any relative neighborhood of d_ there is a vector such that 
3 • < 0. If _z were an element of S, then there would be an 

W 

element r_ m G(Y ) such that ^ = ot£ for some a > 0. However, we 
would then have 


B • z = a 0 ■ r > 0 

which is a contradiction. Therefore z_ is not an element of S and 
is not in the relative interior of S. This completes the proof of 
the proposition. 
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Appendix A 


A Program for Constructing the Natural Cubic Spline Interpolant 


to Given Data. 
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00001 

00002C 

00003C 

00004C 

00005 

00006 
00007 
00008C 
00009C 

0001 oc 
ooouc 
0001 ?c 

00013C 
0001 4C 
00015C 
000 16C 
000 17C 
0001 8C 
00019C 
00020C 
00021C 
00022C 

00023 

00024 

00025 
00026C 
00027C 
00028C 
00029C 

00030 

00031 

00032 C 

00033 100 
00034C 
00035C 
00036C 
00037C 
00033C 
00039C 
00040C 

00041 

00042 

00043 

00044 

00045 

00046 

00047 
00043 200 

00049 

00050 


PROGRAM UNCON ( INPUT -OUTPUT , 1 APL5= INPUT , T APE6=0UTPLIT ) 

WE FORM THE NATURAL CUBIC SPLINE I NTERPOLANT . 

INTEGER N • M . I 

REAL T(50) ,F (50) , IK501 ,X (50), A(50), PP( 4,50) 

REAL AA ( 50 ) , PB < 50 ) , CC < 50 ) 

THE ARRAYS (T) ANB (F) - EACH OF SIZE M, THE NUMBER 
OF DATA . POINTS - CONTAIN THE COMPONENTS OF THE BATA. 
THE DATA FILE IS OF THE FOLLOWING FORM 

M 

T < 1 ) , F (I > 

T< 2) . F(2) 


♦ 

T(M) ,F(M) 

WHERE WE ASSUME (T) HAS STRICTLY INCREASING COMPONENTS. 
READ(3,*) M 

READ(3,*> <T<I),F<I>, 1=1, M) 

N= M-2 

THE ARRAY (D) CONSISTS OF THE SCALED 
SECOND DIVIDED DIFFERENCES. 

DO 100 1=1, N 

D ( I ) = < F( I+2)-F< 1 + 1 ) )/( T( I+2)-T< 1 + 1 ) ) 

- ( F( 1 + 1 )-F( I ) )/( T ( 1 + 1 )-T ( I ) ) 

CONTINUE 


THE SECOND DERIVATIVE OF THE NATURAL CUBIC SPLINE 
INTERPOLANT IS A LINEAR COMBINATION OF LINEAR B-SPLINES. 
WE CALCULATE THE COEFFICIENTS. 


AA < 1 ) = 0.0 

BB < 1 ) = (T<3)-T ( 1 ) )/3 ,0 

CC\1)= (T(3>-T<2) )/6.0 

DO 200 1=2, N-l 

AA< I )= ( T< T + l ) -T ( I > )/6. 0 

BB( I )= (T(I+2)-TU))/3.0 

CC ( I ) = ( T ( I F2) -T ( IF1 ) ) /6 . 0 

CONTINUE 

AA(N)= <T(N+1 >-T(N))/6.0 
BB<N)= (T(N+2)-T(N) 1/3,0 
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00051 

CC(N)= 0.0 

00052 

CALL TRIIi<AA,BB,CC,Ii,N> 

00053C 


00054C 


00055C 


00056 

A<1) = 0,0 

00057 

A(M) = 0.0 

00056) 

do 300 1 = 2 , tm 

00059 

A ( I ) - IKI-l) 

00060 300 

CONTINUE 

00061C 


00062C 


00063C 

NOW WE COMPUTE THE NUMBERS PP(J,I) - THE VALUE 

00064C 

OF THE (J-l)ST DERIVATIVE OF THE NATURAL CUBIC 

00065C 

SPLINE INTERPOLANT EVALUATED AT T(I). 

00066C 


00067C 


00068 

DO 400 K=1 ,N+1 

00069 

DF= F (K+l ) -F(N ) 

00070 

DT= T (K+l )-T(K) 

00071 

DA= A(K+1)-A(K) 

00072 

F'P < 4 , K ) = DA/DT 

00073 

F'P(3,K) = A(K ) 

00074 

F’P ( 2 , K ) = DF/DT - <A(K>/2, + DA/6.)*DT 

00075 

PF'( 1 ,K) = F(K) 

00076 400 

CONTINUE 

00077 

F'P ( 4 , M ) = 0.0 

00078 

F'F'(3,M> = 0,0 

00079 

F'P(2,M) = 0.0 

00080 

F'P < 1 » M ) = F(M) 

00081 C 


00082C 


00083C 


00084 

DO 500 K=1,M 

00085 

WRITEC6.450) K ,T(K) , (PP( I ,K> , 1=1,4) 

00086 450 

FORMAT (5X, 15 ,5F14 .6) 

00087 500 

CONTINUE 

0006)80 


00089C 

WE CREATE A DATA FILE FOR PLOTTING THE (JDER)-TH 

00090C 

DERIVATIVE OF THE NATURAL CUBIC SPLINE INTERPOLANT 

0009 1C 

BY EVALUATING IT AT (MM) EQUALLY SPACED POINTS, 

000920 

INCLUDING THE ENDPOINTS. WE ASSUME THAT (JDER) 

00093C 

HAS VALUE 0, 1, 2, OR 3. 

00094C 


00095 

JDER= 0 

00096 

MM= 201 

00097 

CALL DATAFLU , F’P, M, MM, JDER) 

00098C 


00099 

STOP 

00100 

END 
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00001 

00002C 

00003C 

00004C 

00005C 

00006C 

00007 

00008 

00009 

00010 
00011 
00012 
000] 3 
0001 4 
00015C 
0003 6C 
00017C 
00018 

00019 

00020 
00021 200 

00022 300 

00023 

00024 
00025C 
0002 <LC 
00027C 
00028C 

00029 

00030 

00031 

00032 

00033 

00034 

00035 400 

00036 

00037 450 

00038 500 

00039 

00040 


SURROUTINF DATAFL<TX,PP,LI.HM,JDER) 

WE CREATE A DATA FILE FOR PLOTTING THE (JDER)-TH 
DERIVATIVE OF THE PIECEWISE CUBIC POLYNOMIAL. WE 
ASSUME ( JDER) HAS VALUE 0, 1, 2, OR 3, 

INTEGER LltMMtJDER 
REAL TX( 100) ,F‘F'(4, 100) 

LEFT= 1 
MMONE= MM - 1 
WRITE 1 ! 4,^) MM 

XE= < TX(LI >-TX( 1 ) ) /FLOAT (MMONE ) 

DO 500 IF-1,MH 

XT= TX( 1 > + XEfcFLOAT ( IP-1 ) 

WE FIND THE INTERVAL IN WHICH THE POINT (XT) LIES. 

IF ( LEFT , NE . LI ) THEN 
DO 200 IS-LEFT .LI-1 
IF ( XT «LT , TX< IS+1 ) ) GO TO 300 
CONTINUE 
CONTINUE 
END IF 
LEFT= IS 

WE NOW COMPUTE THE VALUE OF THE POLYNOMIAL AT 
THE POINT (XT) BY USING MESTED MULTIPLICATION. 

H= XT - TX(LEFT) 

FAC= 4.0 - FLOAT (JDER) 

YT= 0.0 

DO 400 M=4 f JDER+ 1 ,-l 

YT= (YT/FACXH + PP(M,LEFT) 

FAC= FAC - 1.0 
CONTINUE 

WRITE(4,450) XT,YT 
FORMAT (F8.4.E18.9) 

CONTINUE 

RETURN 

END 



00001 

51 IPR0U1 INF. 1 Fl J n ( SUIi , DI AG , SUP , P , N ) 

00002 

INTEGER N, I 

00003 

REAL P iN ) , D I AG < N ) , SUB ( N ) , SUP ( N ) 

00004 

IF (N.LE.l) THEN 

00005 

B( 1 )= B(1)/HIAG(1) 

00006 

RETURN 

00007 

END IF 

00008 

DO 111 1=2, N 

00009 

SUP(I)= SUP( I )./DIAG( 1-1 ) 

00010 

DIAG( I )= DIAG(I) - SUB < I ) #SIJP ( I - 1 ) 

0001 1 

B(I)= B<I) - SUB(I)*B(I-1) 

00012 111 

CONTINUE 

00013 

B ( N ) = B(N)/DIAG<N) 

00014 

DO 222 I=N-1 ,1,-1 

00015 

B( I )= <B<I)-SUP(J)*B(I+1 >)/BIAG(I> 

00016 222 

CONTINUE 

00017 

RETURN 

00018 

END 
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Appendix B 


A Program for Constructing the Shape-Preserving Cubic 
Spline Interpolant to Given Data 
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O000J 

00002C 

00003C 

00004C 

oooor.c 

00006C 
00007C 
00008C 
00009C 
OOO IOC 

ooonc 

00012C 
0001 3C 
00014C 
00015C 
00016 
00017 
00013 
00019 
00020C 
00021C 
00022C 
00023C 
00024C 
00025C 
0002AC 
00027C 
00023C 
00029C 
00030C 
0003 1C 
000 32C 
0003 3C 
00034C 
000? fC 
CO 03 iC 
00037 
C003c 
00039 
00040C 
0004 1C 
000422 
00043C 
00044C 
00045C 
00046 


PROGRAM MAIN - INPUT ,011 1 PUT , TAPF5- INPUT , 1 APE6=0UTPUT ) 

WE COMPUTE A SHAPE -F RESERVING INTERF'OLANT 
TO GIVEN DATA. 


NOTE ON THE SIZE OF THE ARRAYS*. 

THE ARRAYS < T > , (F>, AND (A) MUST BE OF LENGTH 
AT LEAST M, THE NUMBER OF BATA POINTS. THE 
ARRAY (TX) AND THE SECOND COMPONENT OF THE 
ARRAY <F'P'> SHOULD BE OF LENGTH 2M. THE ARRAYS 
(X), (Y>. AND (D) MUST BE OF LENGTH AT LEAST M-2. 
THE ARRAT (ID) MUST BE OF LENGTH AT LEAST M-l. 


REAL T(50),F(50) ,X(50),Y(50) r A(50) 
REAL TX ( 100 ) »PP ( 4 . 100) f TL , TF, , cPS 
INTEGER M,N,ITMAX,1,J,IFLAG,MM 
COMMON D ( 50 ) , I P ( 50 ) 


THE ARRAYS (T) AND (F) - EACH OF SIZE M, THE NUMBER 
OF DATA POINTS - CONTAIN THE COMPONENTS OF THE DATA. 
THE DATA FILE IS OF THE FOLLOWING FORM 


M 

T ( 1 > . F ( 1 ) 

T ( 2 ) r F ( 2 ) 


(M/ ,F(M) 


WHCF;L Jf ASSUME (T) HAS STRICTLY INCREASING COMPONENTS. 
sLAD- '.4> M 

RE AB (3.*) 0(I),HI), ’ = 1,.1> 

N- M-2 


•:PS‘ f! .. l-OSITU'E NUMBER USED TO TEST FOR 
CONVEX F'.LL i,( NEWTON'S METHOD - SUBROUTINE (ZERO). 
( ITMAX • IS THE MAXIMUM NUMBER OF ITERATIONS 
WHICH UI PERMIT FOR NEWTONS METHOD TO CONVFRGE. 


1 , ')£- 


0004/ 1 7 nn - 2ti 

0004 SC 

C004"C THL u *1 X IS THE KNOT SEQUENCE (T) WITH THE 

00050C ENDPOINTS 1l r,ND TR DELETED. 
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00051C 

00052 

TL= T ( 1 ) 


00053 

TR= T(M> 


00054 

DO 120 1=1, N 


00055 

X(I) = Kill) 


00056 120 

CONTINUE 


00057C 



00058C 

THE ARRAY (ID CONSISTS OF THE SCALED 


00059C 

SECOND DIVIDED DIFFERENCES. 


00060C 
00061 C 

IT IS IMPORTANT THAT WE IDENTIFY DIVIDED DIFFERENCES 


00062C 

WHICH ARE ZERO. THIS MEANS THAT WE MUST COMPARE TWO 


00063C 

FLOATING-POINT NUMDERS. TO DO THIS WE ASSUME IKK) IS 


00064C 

ZERO IF IKK) IS SMALL. 


00065C 

00066 

XEPS= 1.0 


00067 

DO 130 J=l,20 


00068 

XEPS= XEPS/10 . 


00069 

Z= 1.0 + XEPS 


00070 

IF ( Z .EG. 1.0 ) GO TO 135 


00071 

YEPS= XEPS 


00072 130 

CONTINUE 


00073 135 

CONTINUE 


00074 

YEPS= YEF'S*1000. 


00075C 

00076 

DO 140 K=1,N 


00077 

IK K 1 = ( F (K+2) -F (K+l ) )/( T (K+2)-T (K+l ) ) 


00078 

C - ( F (K+l )-F(K) )/( T (K+l )-T(K) ) 


00079 

IF ( ABS(IKK) ) .LE. YEPS ) IKK)= 0.0 


00080 140 

CONTINUE 


00 08 1C 
00082C 
00083C 

THE INITIAL GUESS (Y) FOR NEWTON'S METHOD 


00034C 

WILL YIELD THE SECOND DERIVATIVE OF THE 


00085C 

NATURAL SPLINE SOLUTION, EXCEPT POSSIBLY 


00086C 

WHEN IK K ) = 0.0 FOR SOME K. 


00087C 

00088 

DO 145 K=) , N 


00089C 

00090C 

00091 

IF ( IKK) .GT, 0,0 ) THEN 


00092 

Y < K > = 1.0 


00093 

ELSE 


00094 

Y(K)= -1.0 

. 

00095 

END IF 


00096C 
00097C 
00098 145 

CONTINUE 


0009?C 

00100C 

00101 

WRITE (6, 150) 
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00102 150 

00103 

00104 160 

00105 

00106 170 
00107C 
00108C 
00109C 
00110C 
0011 1C 
00112C 
00113C 
00114C 
001 15C 
00116C 
001 17C 
00113C 
001 19C 
00120C 
00121 C 
00122C 

00123 

00124 


FORMAT//, ' DATA VALUES ',/) 
WRITE < 6, 160) (IK I ) , 1 = 1 ,N) 
FORMAT (5X,4E12.6) 
WKITE/6,170) 

FORMAT'//) 


I IK K ) = 1 INDICATES THAT THE INTERPOLATING FUNCTION 
IS CONSTRAINED TO BE CONVEX ON CT(K),T(K+1)D 
AND, HENCE, ITS SECOND DERIVATIVE IS CONSTRAINED 
TO BE NONNEGATIVE ON THIS INTERVAL . 

I D ( K ) = -1 INDICATES THAT THE INTERPOLATING FUNCTION 
IS CONSTRAINED TO BE CONCAVE ON CT<K) ,T(K+1 ) 3 
AND, HENCE, ITS SECOND DERIVATIVE IS CONSTRAINED 
TO BE NONPOSITIVE ON THIS INTERVAL. 

I D < K ) = 0 INDICATES THAT THE INTERPOLATING FUNCTION 
IS UNCONSTRAINED ON IT ( K ) , T (K+l ) 1 . 


DO 180 1=1, N-l 
IIKI41) = 0 


00125 

00126 

00127 180 

00128 

00129 

00130 

00131 
0013? 
00133C 

00134 

00135 
00126 
00137 
OC138 
00139C 
001 40C 
00 1 4 1 C 
001 42C 
00143C 

00144 

00145 

00146 185 
00147C 
00140C 
00149C 
001 50C 
00151C 
00152 


IF (D(I).GE.O.O .AND. IK 141 ) .GE.O.O) IIK 1+1 )= 1 
IF (IKI).LE.O.O .AND. IK 14-1 ) .LE. 0.0) IIKI+1)= -1 
CONTINUE 

IF < D(l) .GE. 0.0 ) THEN 
UK 1 )= 1 

ELSE 

IIK 1 )= -1 
END IF 

IF ( D(N) .GE. 0.0 ) THEN 
IIKN+l )= 1 

ELSE 

IIKN+l ) = -1 
END IF 

IF A NONZERO DATA VAlUF IK I ) LIES BETWEEN TWO 
ZERO DATA VALUES IKI-D AND IKI+1), THEN IK I) 

IS TAKEN TO BE ZERO FOR COMPUTATIONAL PURPOSES. 


DO 185 1=2, N-l 

IF < DU-1). EQ. 0.0 .AND. D( 14-1 ) .EG. 0.0 ) IKI)= 0.0 
CONTINUE 

SUBROUTINE /ZERO) CALCULATES THE PIECEWISE 
LINEAR SECOND DERIVATIVE OF THE SHAPE- 
PRESERVING INTERPOLANT, 

CALL ZEF:0( Y,X,N, ITMAX , EPS , IFLAG,TL ,TR) 
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00153C 

00164 

00165 

001 56 

00157 
00153 190 
00159C 
00160 
00161 200 
00162 

00163 210 

00164 

00165 220 
00166C 
001 67C 
00168C 
001 69C 
00170 
00171C 

00172 

00173 230 

00174 

00175 

00176 240 

00177 250 

00178 

00179 260 

00180 
00181 270 
00182C 
00183C 
00184C 
001 BSC 
00186C 
00187C 
00183C 
00189C 

00190 

00191 

00192 
001 93C 

00194 

00195 


6 ( 1 )= 0.0 
A>.M) = 0.0 
DO 190 1=2, N-M 
A ( I ) = Y ( [-1 > 

CONI 3 NUE 

WRITE(6,200) 

FORMAT!/, ' PIECEWISE LINEAR 2ND DERIVATIVE ',/> 

WRITE (6, 210) U,T(I),A(I), 1 = 1, M) 

FORMAT (5X,I5, ' ( ',F14.6,' , *,F14.6,' )') 

WRITE (6, 220) 

FORMAT!//) 

SUBROUTINE (POLY) INTEGRATES THE RESULT 
FROM SUBROUTINE (ZERO). 

CALL POLY ( A,T ,F’F',M,F,LI ,TX) 

WRITE (6,230) 

FORMAT < / , ' KNOTS AND COEFFICIENTS OF PIECEWISE CUBIC',/) 
DO 250 1=1, LI 

WRITE(6,240) I ,TX( I ) , ( PP ( J , I ) , J=l,4) 

FORMAT ( 5X , 15, 5F14 ,6 ) 

CONTINUE 

WRITE (6, 260) I FLAG 

FORMAT (/,' ERROR CODE = ',15,/) 

WRITE (6, 270) IT MAX 

FORMAT!/,' NUMBER OF ITERATIONS =',I5,) 

SUBROUTINE (DATAFL) IS USED TO CREATE A 
DATA FILE FOR PLOTTING. WE EVALUATE THE 
(JDER)-TH DERIVATIVE OF THE PIECEWISE CUBIC 
POLYNOMIAL AT MM EQUALLY SPACED POINTS, 

INCLUDING THE ENDPOINTS TL AND TR. WE 
ASSUME ( JDER) HAS VALUE 0, 1, 2, OR 3 . 

MM= 201 
JDER- 0 

CALL DA r AFL ( TX , PP , LI , MM , JDER > 

STOP 

END 
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00001 

00002C 

00003C 

00004 

00005 
00004 
00007 

ooooa 

00009C 

0001 oc 

0001 1C 
0001 2C 
000 13C 
00014C 
00015C 
00016C 
00017C 

oooi ac 

0001VC 

00020C 

00021C 

00022C 

00023C 

00024C 

00G25C 

00026C 

00C27C 

0002SC 

00029C 

00030C 

90031C 

O0032C 

G0Q33C 

00034C 

00035L 

00037 100 

v0038 C 

00039 

C0O40C 

0004 1 C 

0004 ?C 

00043C 

00044C 

00045 ': 

00046C 

0004 ; c 

00043C 

u0049C 

ooor.or 


SUBROUT INF ZERO ( A , X , N , I THAX , EFS , I FLAG , TL f TR ) 


1NTEGFR N , I TMAX , K , J t L J , L , IFL AG 
REAL A(N),X(N>,FX<50),AL,XL,AR,XR,DT,DA,T,UI 
REAL SUB<50> , HI AG (50) ,3UF‘(50) ,H<50) ,SUM1,SUM2 
REAL RATIO, CLEFT, GRIGH, EPS, FN0RM1,TL,TR 
COMMON B< 50), ID (50) 

INPUT parameters: 

A... INITIAL ESTIMATE FOR NEWTON'S METHOD. 

X. . .KNOT SEQUENCE WITH THE ENDPOINTS DELETED, 

N. . ,THE SIZE OF I HE ARRAY ( A) J I HE NUMBER OF UNKNOWNS, 

UMAX... MAXIMUM NUMBER OF ITERATIONS FOR NEWTON'S METHOD. 

EPS... PARAMETER USED TO TEST FOR CONVERGENCE. 

TL,TR, , .LEFT- AND RIGHT-ENDPOINTS OF THE 
INTERVAL RESPECTIVELY. 

OUTPUT PARAMETERS*. 

A... THE CALCULATED ZERO IF CONVERGENCE OCCURRED. 

UMAX... NUMBER OF ITERATIONS REQUIRED FOR NEWTON'S 
METHOD TO CONVERGE. 

IFLAG, , , JFLAG= 1*. CONVERGENCE INDICATED BY COMPARING 
THE LI NORMS OF THE ITERATES 
TFLAG= 2J NUMBER OF ITERATIONS EXCEEDED ITMAX . 


PRINT 100 

FORMAT ( ' ITERATION NUnBER AND RESIDUAL’.',/ 

' QUADRATIC CONVERGENCE IS EXPECTED,',/) 

DO 350 L J- 1 > ITMAX 

THE ARRAYS (SUB), (IHAG), AND (SUP) CONTAIN 
THE ELEMENTS OF THE TRIDIAGONAL POSITIVE-DEFINITE 
JACOBIAN MATRIX (J), EVALUATED AT THE VECTOR (A). 

IT SHOULD BE NOTED- THAT THE MATRIX EQUATION SOLVER, 
THE SUBROUTINE (TRIP), DOES NOT TAKE ADVANTAGE OF 
THE SiflMETRY OF (J), HENCE (SUB) AND (SUP) ARE 
DOTH NECESSARY. ALTHOUGH SUB(\)=SUP(K-1 ) , EQUATIONS 
FOR DOTH ARRAYS ARE WRITTEN OUT IN FULL. 

I- D ( K >=0.0 FOR SOME K, THEN THE NUMBER 



0005 1C 

OF UNKNOWNS (AND EQUATIONS) REDUCE. IN ORDER 

00052C 

TO PERMIT THE COMPUTATION OF ONE 

JACOBIAN 

00053C 

MATRIX THE PROGRAM SETS SUB ( K ) =SUP ( K-l ) = 0 » 0 

000548 

00055C 

AND DIAG(K)=1 .0. 


00056 

000578 

DO 125 K=1 r N 


00053 

IF (K.EQ.l ) THEN 


0005? 

AL = 0.0 


00060 

XL= TL 


00061 

ELSE 


00062 

Al = A (K-l 1 


00063 

XL= X(K-l) 


00064 

00065C 

00066C 

END IF 


00067 

IF (K.EQ.N) THEN 


00068 

AR= 0.0 


0006? 

XR= TR 


00070 

ELSE 


00071 

AR= A < K+3 > 


00072 

XR= X ( K + l ) 


00073 

00074C 

00075C 

END IF 


00076 

IF ( AL.GE.0.0 .AND. A(K).GE.0.0 ) 

Jl = 1 

00077 

IF ( AL.LT.0.0 .AND. A(K).GE.0.0 ) 

Jl= 2 

00078 

IF ( AL.GE.0,0 .AND. A(K).LT,0.0 ) 

Jl = 3 

00079 

00080C 

IF ( AL.LE.0.0 .AND. A(K).LE.0.0 ) 

Jl= A 

00081 

IF ( A(K) .GE.0.0 .AND. AR.GE.0.0 ) 

J2= 1 

00082 

IF ( A(K) .LT.0.0 .AND. AR.GE.0.0 ) 

J2= 2 

00083 

IF ( A(K). GE.0.0 .AND. AR, LT.0.0 ) 

J2= 3 

00084 

00085C 

IF ( A ( K ) , LE . 0 , 0 .AND, AR.LE.0.0 ) 

J2= A 

00086 

DT= X(K)-XL 


00087 

00088C 

DA= A(K)-AI_ 


0008? 

00090C 

IF ( I D ( K ) .EQ, 1) THEN 


00091 

00092C 

IF (N.NE.l) THEN 


00093 

00094C 

IF (Jl.EQ.l) THEN 


00095 

SUB ( K ) = DT/6,0 


00096 

00097C 

GLEF T = DT/3,0 


00098 

000998 

ELSE IF (J1.EQ.2) THEN 


00100 

T = XL- ( DT /DA > *Al. 


00101 

W= 0,5!) ( X(K >+T ) 
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00102 

00103 

00104 
001 Of. 
00106C 
00107 
00108C 

00109 

00110 
0011 1 
00112 

00113 

00114 
001 3 5C 
00116 
00117C 
00118 
00119 
00120C 
00121 
00122C 
00123 
00124C 


SUB(K)= (X(tO-T)/6,0 * ( ( (l-xt >/DT)*( (X < Kl-T) /PT > 
C + 4.0*( (U-XL>/[iTU( <X(K)-U)/DT) ) 

GLEFT= <X(M-T>/6.0 * ( < (T-XL)/DT)W2 
C + 4.0*<<<W-XL)/DT>**2> + 1.0 ) 


ELSE IF (J1.EQ.3) THEN 


T- XL- ( DT/DA ) ML 
W= 0.5*( T+XL 1 

BUB ( K ) = (T-XL ) / 6 . 0 * ( 4 ,0*( (W-XL)/DT>*< <X(K)-W VDT) 
C + t (T-XL)/DTU((X(M-T)/DD ) 

GLEFT= U-XD/6.0 * ( 4 . 0* < ( < W- XL ) /DT )**2) 

C + <(T-XL)/DD**2 > 


ELSE IF <Jl,En.4> THEN 


SUB < K 1 - 0.0 
GLEFT= 0.0 


ENIi IF 

ELSF IF (K.EQ.H THEN 


00125 

00126 
00127 


00146C 

00147 

00148 
001 v? 
00150 
ooir.i 


SUB( 1 )= 0.0 
GLEFT= 0,0 

IF (Jl.EQ.il GLEFT= DT/3.0 


00128C 

0012? 

00i30E 

ENIi IF 



00131 
001 32C 

ELSE IF ( ID(K) 

, ED. 0 ) 

THEN 

.0133 

SUP ( K 1 - DT/6.0 



00174 
0013 C 

GLEFT= DT/3,0 



O0!3c. 

ELSE IF ( IIKK) 

.EQ. -1 ) 

THEN 

0013" ; C 
L '138 
G0139C 

IF (K.NE.l) 

THEN 


00140 
001 41 C 

IF <J1. 

EQ.4) (HEN 


0C1 4 2 

SUB (K ) = DT/6.0 



00)43 

GLEFT- 1.1/ 3.0 



00144C 

00145 

ELSE IF 

(J1.EQ.3) 

THEN 


T- XL-(DT/DA)*AL 
W= 0.5*< X(KHT > 

SUB ( K ) - ( X ( K ) -T ) /6 . 0 * ( ( (1 -XL ) /DT ) * ( (X <K ) -T)/D1 ) 
D + 4.0*<(U-XL)/DT>*((X(K)-WVDT) ) 

EL EF1 - », -M-TJ/6.0 * \ \ \ T-XL 1 /DT 1 #*2 

C + 4,0*(<<W-XL)/DT)**2> + 1.0 ) 





00 


ooi r.3c 




00154 


ELSE IF (J1.E0.2) THEN 


00155C 




00156 


T= XL-(HT/DA)»AL 


00157 


W= 0.5*< T+XL ) 


001511 


SUB ( K ■> — (T-XD/6.0 * < 4 ,0*( (U-XL) /IU )* ( ( X (K )-U)/DT ) 


0015? 

C 

+ ( ( r-XL)/DT>*( (X(N)-T)/HT) ) 


00160 


GLEFI = (T-XD/6.0 * ( 4.0*( ( <W-XL)/DT)**2> 


00161 

C 

+ <<T-XL)/DT>**2 ) 


00162C 




00163 


ELSE IF (Jl.EB.l) THEN 


00164C 




00165 


SUB(D = 0.0 


00166 


GLEFT= 0.0 


00167C 




00168 


END IF 


00169C 




00170 


ELSE IF (K.EQ.l) THEN 


00171C 




00172 


SUB( 1 )= 0.0 


00173 


GLEFT= 0.0 


00174 


IF (Jl.Ed.4) GLEFT= IH/3.0 


00175C 




00176 


END IF 


00177C 




00178 


END IF 


00179C 




00180 


IF (K.NE.l) THEN 


00181 


IF v D(K-l) .Ed. 0.0 ) THEN 


00182 


SUB ( K ) = 0.0 


00183 


GLEFT = 0.0 


00184 


END IF 


00185 


END IF 


00186C 




00187 


DT= XR-X(K) 


00188 


DA= AR-AOO 


00189C 




00190 


IF ( ID(K+1 ) .Ed. 1 1 THEN 


00191C 




00192 


IF (K.NE.Nl THEN 


00193C 




00194 


IF ( J2.EQ.1 ) THEN 


00195C 




00196 


SUP(N)= DT/6.0 


00197 


GRIGH= DT/3.0 


00198C 




0019? 


ELSE IF (J2.EQ.2) THEN 


00200C 




00201 


T= X(D-<DT/rtA)*A(K> 


00202 


W= 0.5*< XRFT ) 


00203 


SUF'(K 1= (XR-T 1/6.0 * ( < (T-X(K) )/DT)*( (XR-T)/DT) 
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oo;*04 

00705 
00206 
00207C 
00203 
00209C 
00210 
00211 
00212 
00213 
0021 4 
00215 
0023 6C 
00217 
0021SC 

00219 

00220 
0023 1C 
00222 
00223C 
00224 
00225C 
00276 

00227 

00228 
00229C 
00230 
0O231 C 
00232 
00233C 

00234 

00235 
00236C 
00237 
00238C 
00239 
00240C 
00241 
00242C 

00243 

00244 
00245C 
00246 
00247C 
00243 

00249 

00250 

00251 

00252 

00253 
00254C 


l i 4.o#(<w x < k > ) /r« r ) >r < < xf; ui/nn > 

or i CM- (XR-n/6.0* ( < <xk -i )/ium2 
C + 4.0*( ( (XR-U)/DT)»*2> > 

ELSE IF (J2.EQ.3) THEN 


T= X(K)-<nT/riA)*A(K) 

W- 0,5*< T+X( K ) ) 

SUP(N)= (T-X(M)/6.0 * ( 4.0*((U-X(K))/DT)*((XR-W)/HT) 
C + ( < r-X <K))/DD*<< XR-T ) /D r ) ) 

RRIGH= (T-X(M)/6.0 * ( 1.0 + 4 .0*( < (XR-W)/DT)*#2) 
c + ( (XR-T> /nr m2 > 

ELSE IF (J2.EQ.4l THEN 


SUP ( K > - 0.0 
GRIGH= 0.0 


END IF 

ELSE IF (K.EQ.N) THEN 


SUP(N)= 0.0 
GRIGH= 0.0 

IF (J2.EQ.1) GRIGH= DT/3.0 


END IF 

ELSE IF ( IIKK+1 ) .EQ. 0 ) THEN 

SUP(K>= DT/6.0 
GRIGH= DT/3.0 

ELSE IF < IIKK+1 ) .EG. -1 ) THEN 

IF (K.UE.N) THEN 

IF ( J2.EQ.4) THEN 

SUP ( K ) = DT/6.0 
GF;IGH= DT/3.0 

ELSE IF ( J2.EQ.3) THEN 

T= X(K)-<PT/DA)*A<K) 

W= 0.5*< XF:+T ) 

SUP ( K ) = (XR-T1/6.0 * < ((T-X(K))/DT)*((XR-T)/DT1 
C + 4.0*((U-X<K)>/nT)#(CXR-W)/DT) ) 

GRIGH= (XR-T1/6.0 * ( ( (XR-T) /DTI **2 
C + 4.0*(((XR-U)/DT)**2) > 
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0025b 


ELSE IF (J2.E0.2) THEN 

00256C 

00257 


T= X(K)-(Df/riA>*A<K) 

00258 


W= 0.5* < T+X(K) ) 

00259 


SUP ( N 1 - (T-X(M)/6.0 * ( 4.0#< (W-X(K) )/DT)#< (XR-W) 

00260 

C 

+ ( ( T-X(K) )/DT>*( ( XR-T ) /DT ) ) 

00261 


GRIGH= (l-X(M)/6.0 * < 1.0 + 4 .0*< < < XR-U) /DT >##2 

00262 

C 

+ ( ( XR-T )/DT )**2 ) 

00263C 

00264 


ELSE IF ( J2.E0.1 ) THEN 

00265C 

00266 


SUP < K ) = 0.0 

00267 


GRIGH= 0.0 

00268C 

00269 


END IF 

00270C 

00271 


ELSE IF (K.EQ.N) THEN 

00272C 

00273 


SUP(N)= 0.0 

00274 


GRIGH= 0.0 

00275 


IF < J2.E0.4l GRIGH= DI/3.0 

00276C 

00277 


END IF 

00278C 

00279 


END IF 

00280C 

00201 


IF (K.NE.N) THEN 

00282 


IF ( D ( K + 1 ) .EG. 0.0 ) THEN 

00283 


SUP ( N ) = 0.0 

00284 


GRIGH =0.0 

00285 


END IF 

00286 


END IF 

00287C 

00288 


D I AG ( K ) = GLEFT+GRIGH 

002890 

00290C 

00291 


IF ( IKK) .EG). 0.0 ) THEN 

00292 


DIAG(K)= 1.0 

00293 


SUB ( K ) = 0,0 

00294 


SUP ( K ) = 0.0 

00295 


END IF 

00296C 
00297 125 


CONTINUE 

00298C 

00299 


DO 150 L=1,N 

00300 


H(L)= D(L ) 

00301 150 


CONTINUE 

00302C 

00303C 

00304C 


UE SOLVE THE MATRIX EQUATION JX=H, THE ARRAY <H) 
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00305C BEING IDENTICAL TO THE ARRAY (ID. THF. SOLUTION 

00306C IS RETURNED IN THE ARRAY <H>. 

00307C 

00308C 

0030? TALL TRILKSUB,DIAG f SUP,H,N> 

003 IOC 

00311 SUM1 = 0,0 

00312 DO 200 L=1 »N 

00313 A(L) = H(L) 

00314 SUM1 = SUM1 + ABS(A<L>> 

00315 200 CONTINUE 

0031 AC 

00317C THE FUNCTION EVALUATION SUBROUTINE COMPUT MAY 

0031SC BE DELETED. IN THIS CASE THE FOLLOWING EIGHT 

00319C LINES ARE TO BE DELETED AND THE ARRAY (FX) 

00320C CAN BE TAKEN FROM THE REAL STATEMENT AT THE 

00321 C BEGINNING OF THIS SUBROUTINE. 

00322C 

00323 CALL CDMPUT(A,FX,N,X,TL,TR) 

00324 FN0RM1= 0.0 

00325 DO 250 L=1,N 

00326 FN0RM1= FN0RM1 + FX < L > *FX ( L > 

00327 250 CONTINUE 

00328 FN0RM1 = SORT < FN0RM1 ) 

0032? WRITE ( 6 >300) LJ,FN0RM1 

00330 300 FORMAT (15? El 5. 6) 

00331C 

00332C 

00333 IF (LJ.NE.l) THEN 

00334 RAT 10= ABS(SUM1-SUM2> 

00335 AB= EPS&SUM2 

00336 IFLAG= 1 

00337 IF (RATIO .LE. AB) GO TO 400 

00338 END IF 

00339 SUM2= SLIM1 

00340 350 CONTINUE 

00341 IFLAG= 2 

00342 400 CONTINUE 

00343 ITMAX= LJ 

00344 RETURN 

00345 END 
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00001 

00002C 

00003C 

00004C 

OOOOf.C 

00006C 

0000 ? 

00008 

00009 

00010 
00011 
00012C 
00013 
00014C 

00015 

00016 

00017 

00018 

00019 

00020 
00021 
00022C 

00023 

00024 

00025 

00026 

00027 

00028 
00029 
00030C 
00031C 
00032C 

00033 

00034 

00035 

00036 
00037C 

00038 

00039 

00040 

00041 
00042C 
00043 
00044C 
00045 
00046C 
00047 
00048C 
00049 
00050C 


SUBROUTINE COMPUK A,FX,N,X,TL,TR) 

SUBROUTINE (COHPUT), THE FUNCTION EMULATING 
SUBROUTINE , IS OPTIONAL, 


REAL A(N),FX(N),F1,AL0,AHI,TL0»THI 
REAL GLEF,GF:IG,TS,X(N) 

INTEGER N,K,J1 ,J2 
COMMON IK 50) , I IK 50) 

DO 100 K=1 f N 

IF ( IKK) ,NE , 0,0 ) THEN 

IF (K.EQ.l) THEN 

ALO= 0.0 
T L0= TL 

ELSE 

ALO- A(K-l) 

TLO= X(K-l) 

END IF 

IF (K.EG.N) THEN 

AHI= 0.0 
THI= TR 

ELSE 

AHI = A(K+1 ) 

THI= X(K+1 ) 

END IF 


IF (ALO. GE, 0.0 .AND, A(K).GE.O.O) Jl= 1 
IF (ALO. LT. 0.0 .AND, A(K).GE.O.O) Jl= 2 
IF (ALO. GE. 0.0 .AND. A(K).LT.O.O) Jl= 3 
IF (ALO. LT. 0,0 .AND. A(K).LT.O.O) Jl= 4 

IF (A(K).GE.O.O .AND. AHI.GE.O.O) J2= 1 
IF (A(K).LT.O.O .AND, AHI.GE.O.O) J2= 2 
IF (A(K).GE.O.O .AND, AHI.LT.O.O) J2= 3 
IF (A(K).LT.O.O .AND. AHI.LT.O.O) J2= 4 

DT= X (K)-TLO 

IF ( I IKK) .EQ. 1 ) THEN 
IF (Jl.EQ.l) THEN 
GLEF = DT#( 2,Q*A(K) + ALO )/6.0 



00051 

00052C 

00053 

00054 

00055 
00056C 
00057 
00058C 

00059 

00060 
00061 C 
00062 
00063C 
00064 
00065C 
00066 
00067C 
00068 
00069C 
00070 
00071C 
00072 
00073C 
00074 
00075C 
00076 
00077C 
00078 
00079C 
00080 
00081 
00082 
00083C 
00084 
00085C 
00086 
00087 
00089C 
00089 
00090C 
00091 
00092C 
00093 
00094C 
00095 
00096C 
00097 
00098C 
00099 
00100C 
00103 


ELSE IF (J1.EQ.2) THEN 

TS= TLO - ALO*Dl/< ACM-ALO ) 

Fl= ( (TS+X(K) >*0.5 - TLO >/Df 

GLEF= <X(K)-TS)*A<K)*< 2.0*F1 + 1.0 )/6.0 

ELSE IF (J1.EQ.3) THEN 

TS= TLO - ALO*DT/< A(K)-ALO ) 

GLEF= ( <TS-TL0>**2 )*AL0/(6.0*DT) 

ELSE IF (J1.EQ.4) THEN 


GLEF= 0.0 


END IF 

ELSE IF ( irt(K) .EG. 0 ) THEN 

GLEF= DT*< 2.0*A<K> + ALO )/6.0 

ELSE IF ( ID(K) .EG. -1 ) THEN 

IF (J1.EG.4) THEN 

GLEF= DT*< 2.0)fcA<K> + ALO )/6.0 

ELSE IF (J1.EG.3) THEN 

IS= TLO - ALO*DT/< A<K)-ALO ) 

Fl= < (TSfX(K)>*0.5 - TLO >/DT 

GLEF= <X<K)-TS)*A<K)*( 2.0*F1 + 1.0 )/6.0 

ELSE IF (J1.EG.2) THEN 

TS= TLO - ALO*DT/( A(K)-ALO ) 

GLEF= ( <TS-TL0)**2 >*AL0/(6.0*DT> 

ELSE IF (Jl.EG.l) THEN 

GLEF= 0.0 

END IF 

END IF 

DT= THI-X(K) 

IF ( IIKK+1) .EG. 1 ) THEN 
IF (J2.EG.1) THEN 
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00102C 

00103 GRIG= DT*( 2.0*A(K) + AHI )/6.0 

00104C 

00105 ELSE IF (J2.EG.2) THEN 

00106C 

00107 IS= X ( K ) - A(K)*DT/( AHI-A(K) ) 

00108 GRIG= ( (THI-TS>**2 )*AHI/<6.0*DT) 
00109C 

00130 ELSE IF (J2.EG.3) THEN 

00111C 

00112 TS= X(K) - A(K)*DT/( AHI-A(K) ) 

00113 Fl= ( THI-0.5*( TSFX(K) ) >/DT 

00114 GRIG= ( TS-X ( K > )*A(K>*< 1.0 + 2.0*F1 

001 15C 

00116 ELSE IF ( J2.EQ.4) THEN 

00117C 

00118 GRIG= 0.0 

00119C 

00120 END IF 

00121C 

00122 ELSE IF ( IIKK+1) .EG. 0 ) THEN 

00123C 

00124 GRIG= DT*< 2.0*A(K) + AHI 1/6.0 

00125C 

00126 ELSE IF ( IDCK+1) .EG. -1 ) THEN 

00127C 

00128 IF ( J2.EQ.4) THEN 

00129C 

00130 GRIG= DT*< 2.0*A(K) + AHI )/6,0 

00131C 

00132 ELSE IF (J2.EG.3) THEN 

00133C 

00134 1S= X (K) - A(K)*DT/( AHI-A(K) ) 

00135 GRIG= ( <THI-TS>**2 ) *AHI/ ( 6 . O^DT ) 
00136C 

00137 ELSE IF (J2.EG.2) THEN 

00138C 

00139 TS= X < K ) - A(K)*DT/< AHI-A(K) ) 

00140 Fl= ( THI-0.5*( TSFX(K) ) )/DT 

00141 GRIG= ( TS-X(K) )*A(K>*< 1.0 + 2.0*F1 

00142C 

00143 ELSE IF (J2.EG.1) THEN 

00144C 

00145 GRIG= 0.0 

00146C 

00147 END IF 

00148C 

00149 END IF 

00150C 

00151 C 

00152 IF (K.NF..3) THEN 


)/ 6,0 


)/ 6.0 
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00153 

00154 
00155C 

00156 

00157 

00158 
00159C 
00160 
00161(3 
00162 

00163 

00164 
00165(3 
00166 100 

00167 

00168 


IF ( IKk-D .EG. 0.0 ) GLEF= 0.0 
END IF 


IF (K.NE.N) THEN 

IF ( IKM-1) ,EQ. 0.0 ) GRIG= 0.0 
END IF 

FX(K)= GLEF + GRIG - IKK) 


ELSE IF ( IKK) 
END IF 


.EQ. 0.0 ) THEN 

FX(K)= 0.0 


CONTINUE 

RETURN 

END 
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00001 

00002C 

00003C 

00004C 

00Q05C 

00006C 

00007C 

00008C 

00009C 

00010C 

0001 1C 

00012C 

00013C 

00014C 

00015C 

Q0016C 

00017C 

00013C 

00019C 

00020C 

0002 1C 

00022C 

00023C 

00024C 

00025C 

00026C 

00027C 

00028 

00029 

00030 

00031 

00032 

00033 

00034 

00035 

00036 

00037 
00038C 

00039 

00040 

00041 

00042 

00043 

00044 

00045 

00046 

00047 
00048C 
00049C 
00050C 


SUBROUTINE F'OL'i (A,T,PF\M,F,LI ,1X1 

SUBROUTINE POLY INTEGRATES BACK TWICE THE 
POSITIVE PART OF THE PIECEWISE LINEAR SECOND 
DERIVATIVE WHERE THE DATA SUGGESTS THAT THE 
INTERPOLATING CURVE SHOULD BE CONVEX, THE 
NFGATIVE PART OF THE PIECEWISE LINEAR SECOND 
DERIVATIVE WHERE THE DATA SUGGESTS THAT THE 
INTERPOLATING CURVE SHOULD BE CONCAVE, AND 
THE REMAINING PORTION OF THE PIECEWISE LINEAR 
SECOND DERIVATIVE ON THE TRANSITION INTERVALS, 

THE INTEGRATION YIELDS A PIECEWISE CUBIC 
POLYNOMIAL WITH KNOTS GIVEN BY THE SEQUENCE 
(TX). THIS CUBIC POLYNOMIAL INTERPOLATES THE 
DATA AND ITS COEFFICIENTS ARE DENOTED BY THE 
NUMBERS F’P ( J , I ) - THE VALUE OF THE (J-l)ST 
DERIVATIVE OF THE FUNCTION EVALUATED AT TX(I), 
FOR X SUCH THAT TX( I ) .GE .X .LT , TX < 1 + 1 ) THE VALUE 
OF THE CUBIC POLYNOMIAL IS 

F’P < 1 , 1 ) 

+ PP (2,1) * < X-TX(I) ) 

+ (1/2)PF<3,I) # ( X-TX(I) 1**2 

+ (1/6)PP(4,I) * < X-TX(I) >**3 


INTEGER M, J,L,LI 

REAL A( 50) , T ( 50) , F'FK 4 , 100) , F< 50) , TX ( 1 00 ) ,TAU 
REAL DF ,DT ,DA,C,E 
COMMON IK 50) ,ID(50) 

LI= 1 
MN1= M-l 
DO 100 L=1 ,MN1 
DF= F ( L+l )-F (L ) 

DT= T < L+l ) -T < L ) 

DA= A (L+l ) -A(L ) 

JP= 0 

IF (L.EQ.l) THEN 

IF < IK 1 ) ,EQ, 0.0 ) JP= 1 
ELSE IF (L.EQ.MN1) THEN 

IF ( IKM-2) ,EQ. 0.0 ) JP= 1 

ELSE 

C- IKL-1)*IKL) 

IF ( C .EQ, 0.0 ) JF’= 1 

END IF 
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00051 IF (JP.E0.1) THEN 

00052C 

00053 PF'< 4 »LI )= 0.0 

00054 PP(3,LI>= 0.0 

00055 F‘P ( 2 f L I ) = DF/DT 

00056 PP(1,LI)= F(L ) 

00057 TX(LI )= T<Ll 

00058 LI= LI+1 
00059C 

00060 ELSE IF (JP.EQ.O) THEN 

00061 C 

00062 IF (A(L).GE.O.O .AND. A(L+1 ) .GE.O.O) J= 1 

00063 IF (A(L).LT.O.O .AND. A(L+1 ) .GT.O.O) J= 2 

00064 IF (AiL). GT.O.O .AND, A(L+1 > .LT.O.O) J= 3 

00065 IF <A(L).LE.0.0 .AND. A( L+l > .LE ,0. 0) J= 4 

00066C 

00067 IF ( IIKL) .EG. 1) THEN 

00068C 

00069 IF < J , EG . 1 ) THEN 

00070C 

00071 C= DF/DT - (DA/6.0 + A(L)/2.0)*DT 

00072 PP(4,LI)= IiA/DF 

00073 F'P(3rLI )= A(L) 

00074 PP( 2 >LI ) = C 

00075 PP(1»LI)= F(L) 

00076 TX(LI )= T (L) 

00077 LI= LI+1 
00078C 

00079 ELSE IF (J.EG.2) THEN 

00080C 

00081 TAU= T(L) - A(L)*DT/DA 

00032 C= DF/DT - ( A (L+l )**3)*DT/(6,0*DA*DA> 

00083 F'P ( 4 t L I ) = 0.0 

00084 PP(3,LI)= 0.0 

00085 PP(2»LI )= C 

00086 PP(1,LI>= F(L ) 

00037 PP(4f LI+1 )= DA/DT 

00033 PP(3*LI+1)= 0.0 

00089 PP(2,LI+1)= C 

00090 PP(1,LI+1)= C*(TAU-T(L>> + F(L) 

00091 TX(LI)= T(L) 

00092 TX<LI+1>= TAU 

00093 LI= LI+2 
00094C 

00095 ELSE IF (J.EG.3) THEN 

00096C 

00097 TAU= TIL) - A(L)*DT/DA 

00098 E= F(L) - <A(L)**3>*DT*DT/(6.0*DA*DA> 

00099 C- DF/DT + <A(L)**35*DT/(6.0*DA*DA) 

00100 PP(4.LI)= DA/DT 

00101 F’P (3 f LI )= A(L) 
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00102 

PP(2,LI)= C + A(L>*A(L)*DT*0.5/DA 

00103 

F'F' ( 1 fLI )= F(L) 

00104 

PP(4,LI+1>= 0.0 

00105 

PP < 3 , L I + 1 ) = 0.0 

00106 

PP(2fLI+l)= C 

0010? 

PP( 1 T LI+1 )= C*( TAU-T(L) ) + E 

00108 

TX ( LI ) = T(L) 

00109 

1 X ( LI+1 ) = TALI 

00110 
0011 1C 

LI= LI+2 

00112 

00113C 

ELSE IF (J.EQ.4) THEN 

00114 

PP(4,LI)= 0.0 

00115 

PP(3,LI)= 0.0 

00116 

PP(2,LI>= HF/HT 

00117 

PP(1,LI>= F (L 1 

00118 

T X ( 1 1 ) = T(L) 

00119 

00120C 

LI= LI+1 

00121 

00122C 

END IF 

00123 

00124C 

ELSE IF ( ID ( L ) .EQ. 0 ) THEN 

00125 

C= DF/DT - < DA/6.0 + A(L)/2.0)#DT 

00126 

PP(4,LI)= DA/DT 

00127 

PP<3,LI)= A(L) 

00128 

PP(2,LI)= C 

00129 

PF'< 1 , LI )= F<L) 

00130 

TX(LI )= T (L) 

00131 

00132C 

LI= LI+1 

00133 

00134C 

ELSE IF ( ID(L) .EQ. -1 ) THEN 

00135 

00136C 

IF < J.EQ.4) THEN 

00137 

C= DF/DT - (DA/6.0 + A(L)/2.0)*DT 

00138 

PF* ( 4 » L I ) = DA/DT 

00139 

PF'(3,LI)= A(L) 

00140 

PP(2 f LI )= C 

00141 

F’F'C 1 f LI )= F(L) 

00142 

TX(LI )= T (L ) 

00143 

00144C 

LI= LI+1 

00145 

00146C 

ELSE IF ( J.EQ.3) THEN 

00147 

TAU= T(L) - A(L)*DT/DA 

00148 

C= DF/DT - (A(L + l)*»:3)*DT/(6.0!«iA!»:DA) 

00149 

PF'(4,LI)= 0.0 

00150 

PP(3,LI)= 0.0 

00151 

FP < 2 f L I ) = C 

00152 

PP(1,LI)= F(L ) 



00153 PF'<4,LI + 1> = DA/DT 

00154 PP(3»LI+1>= 0.0 


00155 

00156 

00157 

00158 
0015? 
00160C 
00161 
00162C 

00163 

00164 

00165 

00166 

00167 

00168 

00169 

00170 

00171 

00172 

00173 

00174 

00175 

00176 
001 77C 
00173 
00179C 
00180 
00181 
00182 

00183 

00184 
00135 
00186C 
00187 
00188C 
00189 
00190C 
00191 
00192C 

00193 100 

00194 

00195 

00196 

00197 

00198 
0019? 
00200 


F'P < 2 r L I + 1 ) = C 

PP( 1 f LI+1 )= CXTAU-T <L) ) t F(L ) 
IX<LI )- T(L) 

TX(LI+1 )= TAU 
LI= LH2 


ELSE IF (J.EQ.2) THEN 

TAU= T(L) - A(L iX'DT/HA 

E= F(L) - ( A < L ) **3 ) KDTfcDT / ( 6 . 0*DA¥BA ) 

C= DF/DT + <A(Lm3)*DT/(6.0*DA*DA> 

F'P(4,LI) = DA/DI 

PP(3,LI)= A<L) 

PP < 2 » L I ) = C + A < L ) *A < L ) *DT!fcO . 5/DA 
F'P ( 1 » L I ) = F < L ) 

PP(4,LI+1>= 0.0 
PF‘(3»LI+1 )= 0.0 
PP(2fLI+l )= C 

F'P ( 1 r LI + l )= C*< TAiJ-T(L) ) + E 
TX(LI )= r<n 
IXILI+l )= TAU 
LT= LI+2 


ELSE IF (J.EQ.l) THEN 

F'P <4, LI )= 0.0 
F'F'(3,LI)= 0.0 
F'F'(2»LI ) = DF/DT 
F'P ( 1 , L I ) = F(L 1 
TX(LI )= T<L ) 

LI= 1.1+1 


END IF 

END IF 

END IF 

CONTINUE 
PF'(4,LI)= 0.0 
F'P(3,LI)= 0.0 
PP<2»LI>= 0.0 
F'P ( 1 f L I ) = F<H) 
TX(LI1= KM) 
RETURN 
END 
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